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ABSTRACT 


Numerical investigations on a diffusing S-duct with/without 
vortex generators and a straight duct with vortex generators are 
presented. The investigation consists of solving the full 
three-dimensional unsteady compressible mass averaged 
Navier-Stokes eguations. An implicit finite volume lower-upper 
time marching code (RPLUS3D) has been employed and modified. A 
three-dimensional Baldwin-Lomax turbulence model has been 
modified in conjunction with the flow physics. 

A model for the analysis of vortex generators in a fully 
viscous subsonic internal flow is evaluated. A vortical 
structure for modelling the shed vortex is used as a source term 
in the computation domain. The injected vortex paths in the 
straight duct are compared with the analysis by two kinds of 
prediction models. The flow structure by the vortex generators 
are investigated along the duct. 
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Computed results of the flow in a circular diffusing S-duct 
provide an understanding of the flow structure within a typical 
engine inlet system. These are compared with the experimental 
wall static-pressure, static- and total-pressure field, and 
secondary velocity profiles. Additionally, boundary layer 
thickness, skin friction values, and velocity profiles in wall 
coordinates are presented. In order to investigate the effect of 
vortex generators, various vortex strengths are examined in this 
study. The total-pressure recovery and distortion coefficients 
are obtained at the exit of the S-duct. The numerical results 
clearly depict the interaction between the low velocity flow by 
the flow separation and the injected vortices. 
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CHAPTER 1 


INTRODUCTION AND BACKGROUND 


1.1 S-duct Without Vortex Generators 

The subsonic duct is a feature of the air intake 
propulsion systems for modern aircraft whether the speed of 
the aircraft is subsonic or supersonic. Depending on the 
integration of the engine inlet with the airframe, various 
shaped ducts are employed. The intention of duct design is 
to produce high pressure recovery in order to maintain high 
thrust levels, and low flow distortion consistent with 
stable engine operation. It is common to design ducts to be 
as short as possible because of size and weight 
restrictions. Many aircraft employ curved rectangular, or 
circular shaped ducts with constant or varying cross- 
sectional area in the engine intake systems. For example, 
the Boeing 727, Lockhead Tristar (L-1011) , General dynamics 
F-16, and McDonnell-Douglas F-18, etc., use the S-shaped 
duct in their engine intake systems. Usually, the 
diffusing duct is employed in the inlet of the propulsion 
system of the aircraft in order to decelerate the flow and 
achieve high pressure recovery at the engine compressor. 

The S-shaped duct produces complex cross flow patterns 
and nonuniform velocity profiles at the exit because of its 
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curvature and centerline offset. These deteriorate the 
performance of the engine inlet system. The nonuniform flow 
at the exit results from the expulsion of low velocity 
fluid by a pair of counter-rotating vortices, which are 
produced near the inflection point of the duct and 
stretched toward the exit. 

The experimental results obtained by Bansod and 
Bradshaw ( 1972 ) show the expulsion of low velocity fluid at 
the exit. The authors conducted experiments using three 
different kinds of constant-area S-shaped ducts in 
incompressible flow. The S-shaped ducts were assembled with 
different radii of curvature (R) of the duct centerline. One 
had the same R/D = 2.25 in the first and second half bend. 
Others had R/D = 2.25 or R/D = 3.5 in the first and second 
half bend, respectively. The S-shaped duct with large 
radius in the second half bend was more efficient because 
the thick boundary layer in the second half bend was less 
rapidly deflected. 

McMillan (1982) conducted experiments using a 40° 
curved rectangular diffusing duct. The flow was 
incompressible. The results show a pair of counter-rotating 
vortices at the exit. The secondary velocity profiles show 
that the high velocity fluid at the central portion of the 
channel moves toward the concave wall, driven by 
centrifugal force. Correspondingly, the low velocity fluid 
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in the boundary layers moves toward the convex wall. The 
mean velocity in the diffusing duct is dominated by this 
secondary flow. 

Guo and Seddon(1982) tested a S-shaped rectangular 
duct of constant cross-sectional area with several 
different angles of attack. The results show that the flow 
separation, turbulent intensity, and flow distortion at the 
exit increase with increasing the angle of attack. 

Vakili et al.(1987) tested a diffusing 30°-30° S-duct 
with circular cross section. The duct area ratio between 
inlet and exit was 1.51. The offset of the duct resulting 
from the centerline curvature was 1.34 times the inlet 
diameter. Two straight circular pipes were attached 
upstream and downstream of the S-duct to provide the 
desired boundary layer thickness flow at the inlet of S- 
duct and minimize the exit flow effect. The entrance Mach 
number was 0.6 and the Reynolds number based on the inlet 
diameter was 1.76xl0 6 . The secondary velocity profiles, 
static- and total-pressure contours, and surface static- 
pressure were measured at the several streamwise locations. 
The experimental results show that a pair of counter- 
rotating vortices created by the flow separation cause the 
flow distortion at the exit of the S-duct. 

Jenkins and Loef f ler ( 1991) conducted experiments on a 
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compact diffusing S-duct. The offset of the duct centerline 
was 1.5 times inlet diameter. The duct area ratio between 
inlet and exit was 2.25. The entrance Mach number and 
Reynolds number were 0.34 and 5.75x10 s , respectively. The 
authors measured the secondary velocity profiles, 
streamwise velocity contours, and surface static-pressure. 
The results were similar to the experimental results 
obtained by Vakili et al.(1937). 

Wellborn et al.(1992) conducted experiments on a 
diffusing S-duct, which was larger than, but geometrically 
similar to the duct studied by Vakili et al.(1987). The 
duct inlet Mach number was 0.6 and the Reynolds number 
based on inlet diameter was 2.6xl0 6 . Two straight pipes of 
3.75 times inlet diameter were attached upstream and 
downstream of the S-duct to provide a uniform inflow and 
minimize the exit flow effect. The authors measured the 
surface static-pressure along the streamwise and 
circumferential direction. Streamlines near the wall, 
observed by oil flow visualization, showed the formation of 
the counter-rotating vortices in the flow separation 
region. The results showed that the flow at the exit was 
strongly affected by these vortices, and the mean velocity 
profiles were very similar to the total-pressure field. 

Early numerical work on the curved pipe is shown in 
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Rowe(1970)'s work. The author computed the secondary flow 
on a 4 5 0 -4 5 0 S-shaped and a 180° pipe by the step by step 
application of the Squire and Winter ' s ( 1951) inviscid 
secondary flow theory. The computation based on the 
inviscid theory predicts roughly the flow pattern in a 
curved pipe, if the mean flow does not have large local 
variations associated with the secondary flow. 

Pratrap and Patankar ( 1975) calculated mean velocity 
and secondary flow in a 90° curved constant-area 
rectangular duct for incompressible flow. The authors used 
the fully parabolized Navier-Stokes (PNS) equations and 
partially PNS equations, with a k-£ turbulence model. The 
partially PNS equations for subsonic flow are obtained from 
the full Navier-Stokes (FNS) equations by assuming that the 
streamwise viscous diffusion terms are negligible compared 
to the normal and transverse viscous diffusion terms. The 
fully PNS equations have one more restriction, that the 
pressure in the streamwise momentum equation is assumed to 
vary only in the streamwise direction. More detail 
information about PNS equations is described by Anderson et 
al.(1984). The computational results show that predictions 
using the partially PNS equations are more accurate than 
those using the fully PNS equations. 

Levy et al.(1980) conducted computations in a 
constant-area S-shaped duct using PNS equations. The offset 
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and length of the duct was one and five times the inlet 
diameter, respectively. The inflow Mach number was 0.2. The 
results show that the total-pressure contours at the region 
of near the bottom wall are almost the same shape as the 
streamwise velocity contours. More detailed results of the 
total-pressure contours and secondary velocity profiles 
were obtained by Towne and Anderson ( 1981) . The authors also 
conducted a numerical study with a PNS computer program 
with an algebraic turbulence model. The flow was laminar 
with an entrance Mach number of 0.2 and a Reynolds number 
based on duct diameter of 2000. 

Levy et al. (1983) analyzed a 22.5°-22.5° S-shaped duct 
in laminar and turbulent flow at Reynolds numbers of 790 
and 4.8xl0 4 , respectively, using a PNS computer code with an 
algebraic turbulence model. The streamwise velocity 
contours agreed well with the experimental data. The 
analysis shows that the streamwise velocity in turbulent 
flow is similar to the laminar flow field, but the 
streamwise velocity distortion in the turbulent flow is 
less than that in the laminar flow. 

Vakili et al. (1983, 1984) performed numerical analysis 
and experiments on a 30°-30° non-diffusing S-duct. The 
inlet Mach number was 0.6 and the Reynolds number was 
1.76xl0 6 . The PNS computer code was used to predict the 
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static- and 

total-pressure contours and 

secondary velocity 

profiles . 

The computational 

results 

showed 

that the 

secondary 

velocity profiles 

agreed 

well 

with the 


experimental results. The extent of the flow distortion was 
underestimated due to simplifications made in the pressure 
field calculations. The pressure in the streamwj.se 
direction was used a sum of the pressure obtained from a 
three-dimensional potential flow analysis and one- 
dimensional correction. Harloff et al. (1992a) used the 
three-dimensional FNS equations to analyze the 30°-30° 
nondiffusing S-duct, which had the same geometry and flow 
conditions tested by Vakili et al.(1984). The authors used 
two kinds of grid, H- and O-grid. An H-grid conforms well 
to the rectangular shape. An O-grid, which has a pole 
boundary condition at the center of the grid, conforms well 
to a circular cross-section. The results obtained using the 
O-grid were better than those by H-grid because the H-grid 
had a large amount of grid skewness in the corner region. 
The authors concluded that the computational results were 
in qualitative agreement with the experimental results, and 
more advanced turbulence model and grid refinement could 
improve the agreement with the experimental results. 

Jenkins and Loef f ler (1991) conducted computations on 
a compact diffusing S-duct, and compared their results with 
experimental data. Results were obtained using the Baldwin- 



Lomax and "one-half" equation turbulence model which 
accounts for some of the history effects in computing the 
turbulence length scale. The results showed that the thin 
layer Navier-Stokes equations code provided a reasonably 
good representation of the flow at the exit, but the code 
could not accurately predict the separated flow region. 

Harloff et al. (1992b) conducted a numerical study in 
the diffusing 30°-30° S-duct using the three-dimensional 
FNS equations. The authors used the algebraic and k-£ 
turbulence model. The wall static-pressure distribution and 
total-pressure profiles calculated with the k-£ turbulence 
model were better than those with the algebraic turbulence 
model. However, the computational results showed that both 
turbulence models could not adequately account for strong 
secondary flows with flow separation. 


1.2 S-duct With Vortex Generators 

From the review of the S-shaped duct without vortex 
generators, one sees that the strong secondary flow due to 
adverse pressure gradient may have deteriorating effects on 
the performance of the engine inlet system. To alleviate 
this problem, a vortex generator can be used as a flow 
control device because it can transport energy into the 
boundary layer from the outer flow. The vortex generator 
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has been used mainly for the prevention of separation on 
wings, diffusers, or bends, or at least for decreasing the 
extent of separated region. There are many kinds of vortex 
generators, such as simple plow, shielded plow, triangular 
plow, scoop, twist interchanger, ramp, tapered fin, dome, 
shieled sink, etc.,. Schubauer and Spangenberg (1960) 
experimentally investigated the mixing rate of the 
turbulent boundary layers with many different vortex 
generators in a region of adverse pressure gradient. Most 
vortex generators in use today are small wing sections, 
which are mounted upstream of the problem flow area. The 
vortex generators are inclined to the oncoming flow to 
generate shed vortices. The vortex generators are usually 
sized to local boundary layer height to obtain the best 
interaction between the shed vortex and the boundary layer. 
The vortex generators are usually placed in groups of two 
or more upstream of the problem flow area. Fig. 1.1 shows 
a wing type vortex generator. 

Boundary layer control by vortex generators relies on 
induced mixing between the external or core stream and the 
low energy flow region. The mixing is promoted by 
longitudinally trailing vortices over the duct surface 
adjacent to the edge of the boundary layer. Fluid particles 
with high momentum in the streamwise direction are swept 
along helical paths toward the duct surface to mix with 
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and, to some extent, replace the low momentum boundary 
layer flow. This is a continuous process that provides a 
source of re-energization to counter the natural boundary 
layer growth caused by surface friction, adverse pressure 
gradients, and low energy secondary flow accumulation. 

There are two basic configurations. In one 
configuration, all of the vortex generators are inclined at 
the same angle with respect to the oncoming flow direction, 
as shown in Fig. 1.2(a). These are called co-rotating 
configurations because the shed vortices rotate in the same 
direction. In the other configuration, the vortex 
generators are grouped in pairs inclined in the opposite 
direction to the oncoming flow, as shown in Fig 1.2(b). 
These are termed the counter-rotating configurations 
because the shed vortices in pairs rotate in opposite 
directions to each other. 

What kind of configuration is chosen depends on the 
location of the flow separation for a given geometry. Co- 
rotating vortex generators are very competitive with 
counter-rotating vortices in reducing the flow separation 
if the generators are properly selected and located. This 
type of vortex generator has the following characteristics 
when it is used within the duct. (1) Two induced vortices 
move along the duct surface, (2) the first vortex moves 
away from the duct surface, (3) the other vortex remains 
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close to the duct surface, and (4) the distance between 
them decreases because they counteract each other. 

Counter-rotating vortex generators are very effective 
in reducing the flow separation if the vortex generators 
are placed slightly upstream of the region of separation. 
If the induced vortices are rotating away from each other, 
the induced secondary flow between two vortex generators 
moves toward the center of the duct. The vortices are 
attracted to each other for a short time, and then they 
proceed to march away from the wall. Since the two vortices 
are moving toward the center of the duct, the duct surface 
is not much affected by the induced vortices. If the 
induced vortices are rotating toward each other, the 
induced secondary flow between two vortex generators moves 
toward the duct surface. Two vortices move away from each 
other, but they remain close to the duct surface because 
the induced secondary velocities push each other toward the 
surface. The induced vortex strength is dissipated 
significantly as it moves downstream due to viscous 
diffusion. 

Early studies with vortex generators have focused on 
improving the diffuser performance. Brown et al. (1968) 
conducted experiments with pairs of vane type vortex 
generators in a short diffuser. The results show that high 
pressure recovery and flow uniformity can be achieved by 
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the vortex generators, but the incorrect arrangement of 
vortex generators can lead to significance performance 
losses . 

Vakili et al. (1986) experimentally investigated the 
performance of the vortex generators in a diffusing 30°-30° 
S-shaped duct with circular cross-section. The entrance 
Mach number was 0.6 and the Reynolds number based on the 
diameter was 1.76x10 s . To eliminate the total-pressure 
distortion at the exit and flow separation in the duct, arc 
wing type, rail type and vane type vortex generators were 
installed at the upstream of the separation region. Using 
a flow control device, the flow distortion at the exit was 
significantly improved. The results showed that the flow 
field at the exit depended on the types of vortex 
generators . 

Reichert and Wendt (1992) conducted experiments to 
examine three parameters of vortex generators array, i.e., 
the height of vortex generator, the location of the vortex 
generators array, and the vortex generators spacing. The 
test was performed on the same geometry and flow conditions 
as studied by Wellborn et al.(1992). The Wheeler wishbone 
generators, which produced a pair of counter-rotating 
vortices, were used. The results show that the efficiency 
of vortex generators is much dependent on the parameters of 
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vortex generators. 

A numerical study of a fully viscous subsonic internal 
flow with vortex generators was reported by Kunik(1986). 
The shed vortex is modeled by introducing a vorticity 
source term into a modified form of the PNS equations. That 
vortex model resembles the one proposed by Squire (1965) 
except that it neglects the variation of viscosity in the 
cross plane. Quantitative comparisons with the experimental 
data by Vakili et al.(1986) show that the vortex model can 
predict the global flow field in the S-duct. 

Anderson (1991) conducted the analysis of the flow 
physics associated with vortex injection in the S-shaped 
duct and F/A-13 inlet duct. The author used the PNS 
equations with the algebraic turbulence model. Predicted 
total-pressure profiles were in good agreement with 
experiment results, but the transverse velocities at the 
exit were overestimated. 

The PNS equations were derived from the FNS equations 
using a series expansion technique. These equations can be 
solved using a space-marching technique because the 
streamwise diffusion term in the FNS equations is neglected 
and a pressure in the streamwise momentum equation is 
assumed to vary only in the streamwise direction for 
subsonic flow. Naturally, a substantial reduction in 
computation time and storage is achieved, but the space 
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marching method is not well posed if the streamwise 
pressure gradient is included everywhere in the flow field. 
If the streamwise velocity deficit in the vortex core is 


quickly recovered along the 

duct, 

the 

role of 

the 

streamwise 

diffusion term in 

the 

FNS 

equations 

is 

important. 

In that case, the 

FNS 

equations should 

be 


solved . 



15 



(a) Top view 

U Shed Vortex 



Fig- l.l A typical vortex generator. 




(a) Co-rotating. 



(t>) Counter-rotating. 


Fig. 1.2 Typical vortex generator configurations. 



CHAPTER 2 


GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 


From the reviews of the S-duct without vortex 
generators, we can conclude that computational fluid 
dynamics (CFD) studies have generally used the PNS computer 
code to predict the flow fields in the curved ducts, and 
simple turbulence modelling without modifications cannot 
predict correctly the flow fields which have strong 
secondary flows with flow separation. The PNS solutions 
usually rely on an input inviscid static-pressure field, 
which is generally from an Euler or potential analysis. In 
the present study, the three-dimensional FNS equations with 
a modified algebraic turbulence model are solved to predict 
the flow fields in the diffusing 30°-30° S-duct. The inlet 
Mach number is 0.6 and the Reynolds number based on the 
inlet diameter is 1.76X10 6 . Several aspects of the flow 
fields are examined. The computed static- and total- 
pressure fields, secondary velocity profiles and boundary 
thickness are compared with experimental results obtained 
by Vakili et al. (1986 , 1987) and Wellborn et al.(1992) for 
CFD validation. Additionally, skin friction values and 
velocity profiles in wall coordinates are investigated. 

From the reviews of the S-duct with vortex generators. 
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we know that most of predictions of flow fields with vortex 
generators are conducted using PNS equations. In order to 
apply space-marching method of PNS equations, the vortical 
structures, modelled from the shed vortex, are set up at 
the inlet plane of a computational domain with the 
approximately calculated inlet flow conditions. In contrast 
to the previously published work, a new vortex model is 
developed and it is applied inside the computational domain 
like a source term. The inlet boundary conditions are not 
affected by the added vortical structures. Numerical 
analysis is conducted using the three-dimensional FNS 
equations, with an algebraic turbulence model, because FNS 
equations are able to deal with the streamwise diffusion 
terms, which are important in the region of the shed vortex 
core. In order to confirm the developed vortex model, four- 
different types of vortex generators are examined in a 
straight duct. In the straight duct computations, the inlet 
Mach number is 0.6 and the Reynolds number based on the 
diameter is l.OxlO 6 . The computational results are compared 
with the analytic results obtained by the two prediction 
models. 

In order to investigate the flow structure in the 
diffusing 30° -30° S-duct with vortex generators, the above 
mentioned vortex model is applied inside the computational 
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domain. The three-dimensional FNS equations with a modified 
algebraic turbulence model are solved. The inlet Mach 
number is 0.6 and the Reynolds number based on the inlet 
diameter is 1.76xl0 6 . The interaction between the injected 
vortices and separated flow is investigated. The static- 
and total-pressure fields and secondary velocity profiles 
are compared with the experimental results obtained by 
Vakili et al.(1986). In order to investigate the effects of 
the injected vortices, the computed results are compared 
with those without vortex generators, and the total- 
pressure recovery and distortion coefficients are 
investigated at the exit of S-duct. 


2.1 Governing Equations 

The three-dimensional and compressible Navier-Stokes 
equations in Cartesian coordinates without body forces are 
written in a conservation form as follows: 

dU d(E-E v ) d(F-F v ) d(G-G v ) n (2.1) 

dt dx dy dz 

O is the independent variable, E, F and G are the 


convective flux vectors 
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E v , F v and 


U = 


'p' 

pu 
pv 
p w 
\pe) 


( 2 . 2 ) 


E = 


' pu ^ 
pu 2 +p 
puv 
puw 

k u ( pe+p) j 


(2.3) 


F = 


' pv ' 

puv 
p v 2 *p 
pvw 

,v(p e+p), 


(2.4) 


G = 


pw 

puw 

pvw 

pw 2 +p 


\w(pe+p) ) 


(2.S) 


are the viscous flux vectors: 
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( 2 . 6 ) 


(2.7) 



( 2 . 8 ) 


The first row of the vector Eq. (2.1) corresponds to 
the continuity equation, the second, third and fourth rows 
are the momentum equations, the fifth row is the energy 
equation; e in the energy equation is the summation of 
internal energy and kinetic energy per unit mass. The shear 
and normal stresses can be expressed using Stokes 
hypothesis, i.e., the second viscous factor \=-2n/2 



iu/in + iz ♦ isf\ 

3 H \ dx dy dz I 
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yy 


= 2 V- 


dv _ 2_ 
dy 3 

dw _ _2 
3z 3 


zx xz 


in + 

av + 


3x 

ay 

3z 

Ju + 

av + 

dw 

3x 

ay 

3z 

in + 

3w \ 


3z 

3x / 



(2.9) 


= x = Jir + i"\ 
2y az a y / 


yz 


^jo' ^yx 


= Jin + Jr\ 

^\a y ax/ 


In the energy equation, heat flux g x , g^ and g 2 are 
expressed as; 


a = -*il 

g * 3x 




( 2 . 10 ) 


<7* = 


3z 


To close this system, the state equation with an 
assumption of a perfect gas is employed. 


P = P RT 


(2.11) 


The viscosity, heat conductivity coefficients and 



23 


specific heats of air are evaluated using fifth order 
polynomials in temperature, using properties presented in 
the National Bureau of Standards tables ( 1955) . 

For turbulent flow, it is convenient to use a 
conservation form of the mass-averaged Navier-Stokes 
equation. This form takes all turbulence effects into 
account by adding the eddy viscosity to the equations. 
These equations can be obtained by replacing the molecular 
coefficient of viscosity n with + ji, and also the 
coefficient of thermal conductivity k with k + k t . /x t is the 
eddy viscosity and k t is the turbulent thermal conductivity. 
The turbulent thermal conductivity can be expressed in 
terms of the eddy viscosity and turbulent Prandtl number 
Pr t , i.e., k t = c^/Pr,. In the present study, the turbulent 
Prandtl number is assigned Pr, = 0.9 for air, and the eddy 
viscosity will be discussed in the section on turbulence 
model . 


2.2 Coordinate Transformation 

The computation of flow- fie Ids in and around complex 
shapes such as ducts, engine intakes or aircraft, etc., 
involves computational boundaries that do not coincide with 
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coordinate lines in physical space. For numerical methods, 
the imposition of boundary conditions for such problems has 
required a complicated interpolation of the data on local 
grid lines, and typically a local loss of accuracy in the 
computational solution. Such difficulties motivate the 
introduction of a mapping or transformation from physical 
(x, y, z) space to a generalized curvilinear coordinates 
(£, Tj , f) space. The generalized coordinate domain is 
constructed so that a computational boundary in physical 
space coincides with a coordinate line in generalized 
coordinate space. It makes it possible to solve the 
governing equations on an uniformly spaced computational 
grid. In order to use an uniform grid, consider a general 
transformation of the governing equations. 

£ = l { x, y, z) 

r| = r| ( x, y , z ) (2.12) 

C = C ( x, y , z) 

Using the chain rule of partial differentiation, 
the partial derivatives become 
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± .5 _L + T1 ± + c _i 

ax ^ x ac lx 8n w ac 


_a_ 

8y 


- «. 


ae 


a A , a 
^ an ac 


(2.13) 


_a_ 

3z 


5 JL + n — + C — 

- ac a n z ac 


The Jacobians of the coordinate transformation are as 
follows : 


J = det 


' x t \ V 

n y< 

z i z r i z i } 


(2.14) 


The vector Eq. (2.1) can be written in terms of a 
generalized nonorthogonal curvilinear system (£, 77 , f) 

using the change rule of partial differentiation and the 
Jacobian of the transformation. The resulting equations can 
be written: 

dU d(E-E v ) 8 (/-/„) d(G-G v ) 

at + d( + an + ac 


(2. IS) 
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' pJ(uZ x +vZ y +wl z ) N 
P«7U ( Ul x + V£y+WZ Z ) +PJZ X 
E - p Jv { ui x +v£ y + w£ z ) +pJty (2.17) 

pjV( U^ x + Vt,y+WZ z ) +PJZ Z 
< (pe + p) J(u5/^ y +k^) , 

r Pi7( UT) X + VT| y + vqj > 

P«7U { UT1 X + VT\ y +wr\ z ) +pjT l x 

_P = pJV ( ur\ x + VT\ y + wt| 2 ) +pJT| y (2.18) 

pJV ( L 2 Tj x + VTl y + WTj z ) + p«7r| 2 
k {pe+p) J{ur\ x +vr\ y +wT\ z ) J 

' pj( uC x + v( y +fc< 2 ) ^ 

pJu ( u£ x + v£ y + ^C r ) + P^C X 

G = pcrv ( UC X + ^Cy + W C Z ) + PJCy (2.19) 

P Jv ( u£ x + vC y + vC z ) + pJCz 
v ( pe+p) J( uC x + vC y +wC z ) y 
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o \ 
a vl 



e v2 


£ = 

e v2 

= 


®V4 



, e vS , 

V 


(f v A 

/ 


f V Z 


*v = 


= 





< -^v5 ) 

\ 


'ffvl' 

/ 


9v2 


G v = 

Jv3 

= 


Jv4 ! 



^vsj 

V 


^ XX + X xy^^y + ^ xx^^> z 

+ ^ yy'J^y + ^yz 1 ^"^ z 
^xz^^x + ^yz^y + ^zz^^z 

ue +ve +ve + feii(c — +C — 
Ue v 2 VS v 3 Ve v 4 p r ' C 1 ££ C 2 0 ^ 


c 31 ) 

c 3 *r } 


ac 

( 2 . 20 ) 


V^lx + ^xy 4 ^*! y + ^ JCX 
^xy^^x + ^yy^^y + ^yz^^l z 
^ xz^^x + ^yz'^ly + ^ zz^^ z 

ut„*vf w ,.»t r ..%lc,%.c t %.c,%) 


•V3 


(2.21) 


T JC + t JC + t JC 

xx s x xy > y xx s ; 

T JC + T JC + T JC 

xy ’x yy ^y yz s i 


^xz^^x + ^yz^^y + ^ zz^^ z 

u^^^<c^c s JI^£, 


( 2 . 22 ) 


where 


Cl = jU X 2 + ly 2 +l Z 2 ) 

c 2 = j ( S* 1 )** £ y n y + 5*n*) 

C 3 = J( £xCx + S y C y + $ z C z ) 


C 4 - J { T \ x 2 +T) y 2 +11/ ) 


( 2 . 23 ) 
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C s = J ( T| x C x + T| y C y + T|,C 2 ) 


= *7( C, 


c , 2 - C 2 ) 


Note that all the stress terms in £ v , F» an< * Gv should 
be transformed. For example, the shear stress term would 

be transformed to; 


V = MS 


d u 
y dl 


„ du . r du . c dv , „ dv 

^ y dr\ ^ y ac T 5 x as ar, 


C x 4?) (2.24) 

ac 


2.3 Boundary Conditions 

The numerical solution of any partial differential 
equation requires the application of appropriate number of 
properly posed boundary conditions. The important aspects 
of boundary condition development are that the physical 
definition of the flow problem must be satisfied and the 
numerical algorithm with the developed boundary conditions 
must be stable. The theory of characteristics suggests how 
to decide the conditions required at a boundary. The 
concept of characteristic theory is most easily developed 
for the one-dimensional Euler equations. Extending the 
concepts to three-dimensions, we can obtain three U., U,,, + 
c and U,, - c characteristics in this system, using the fact 
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that the incoming flow at the inlet plane is uniform. At 
the inlet plane, only four pieces of information enter the 
domain along the incoming characteristics and one piece 
leaves along the outgoing characteristics because the flow 
speed in the whole computational domain is subsonic. 
Therefore, four boundary conditions must be specified, and 
one relation has to be extracted from the characteristic 
equation . 

It is not necessary to fix values in terms of the 
actual characteristic variables as long as the alternative 
choice leads to a well posed problem. A particular good 
choice on physical grounds is to specify the stagnation 
enthalpy and the entropy of the incoming flow. For a 
perfect gas, this corresponds to specifying the stagnation 
temperature and pressure. These conditions are same as the 
flow conditions through a duct or nozzle fed from a large 
reservoir in which conditions remain constant. Constant 
stagnation temperature T a and pressure p 0 are specified at 
the inlet plane. 


T = r + — — - 

0 2 c. 


T_ 

To 



( 2 . 25 ) 


Po = P 


( 2 . 26 ) 
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The transverse velocities are assumed to be zero. The 
relation which corresponds to the negative characteristic 
can be derived along the characteristics Eq. (2.27). 


dp 

d$_ 


P c 


d u 
dL 


= 0 


(2.27) 


This becomes 

p - p cu = (p - p cu) ince:ioz (2.28) 

Substituting the inflow boundary condition Eqs. (2.25) 
and (2.28) into Eq. (2.26), we can obtain the inlet 
temperature, static-pressure and axial velocity. The axial 
velocity near the duct edge approaches zero in order to 
satisfy the no-slip condition on the wall during iterations 
because the finite volume method is employed. The density 
at the inlet is calculated from the equation of state. The 
total initial energy at the entry plane can be obtained 
from the calculated values. 

At the exit, one negative characteristic enters 
through the boundary into the computational domain. One 
boundary condition must be specified at this plane. In this 
study, constant static-pressure is specified at the exit 
plane. Physically this condition corresponds to a duct with 
an unobstructed exit into a large constant pressure 
reservoir. Linear extrapolation is adopted for 
evaluating the exit velocity and exit density. The exit 
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temperature is calculated from the equation of state. The 
total internal energy at the exit plane can be obtained 
from the calculated values. 

The no-slip condition is specified on the wall of the 
duct and an adiabatic wall condition is imposed by setting 
the normal derivative of temperature equal to zero. The 
boundary values at the center, which is needed when using 
the O-grid, are evaluated by averaging the surrounding flow 
properties . 


2.4 Turbulence Model 

The Baldwin-Lomax turbulence model (1978) is applied 
along the normal direction from the wall. This model has 
been used extensively for attached or slightly separated 
flows because it leads to low computational time 
requirements and it appears to be comparable to more 
complex turbulence models. The Baldwin-Lomax model is an 
algebraic two-layer eddy viscosity model based on the 
Cebeci-Smith(1974) method with modifications that avoid the 
necessity for finding the edge of the boundary layer. Near 
the wall, the Baldwin-Lomax model uses the well-known 
Prandtl-Van Driest formulation for the turbulent viscosity. 
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<He> 


inner 


= P i 2 


2 i<3 1 


( 2 . 29 ) 


where 


1 


icy [ 1 - exp ( ) ] 

A 


( 2 . 30 ) 


c3| is the magnitude of the local vorticity vector. 



3u _ dv 
By dx 


I dv _ dw 
\ dz By 


dw du \ 2 
dx dz J 


( 2 . 31 ) 


and 


* = = /p^y 


( 2 . 32 ) 


Since the damping constant A is a function of the 

pressure gradient, an empirical equation by Kays and 

+ 

Crawford (1980) is employed for A . 

A * - 

7 .lip* + 1.0 

( 2 . 33 ) 

If p’ > 0 . 0 , b = 2.9 
If p* < 0.0 , b = 4.25 


where 


,(1/2) t _<3/2) 


P 


( 2 . 34 ) 
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In this study, r w is evaluated as the absolute value; the 
value of p + in the computation is less than O(10' 2 ). 

In the outer region, a Clauser formulation with a 
Klebanoff intermittency function is used. 

(M oucar = 0.02688 p F uak ,F klab {y) (2.35) 


where 

^uake ~ ( y i max ^mdx ' 0 • 25 J/' rax Qdif / ^\max ^ 

and Klebanoff intermittency factor is given by 
F kleb (y) = [ 1 + 5.5 ( 0.3 ) 6 ] - 1 

-/max 

where Ymax an< * F max are determined from the equation 

F(y) = y | (3 | [ 1 - exp (--£_)] (2.38) 

The quantity q ^ is the difference between the maximum and 
minimum total velocity in the profile. The parameter is 
the maximum value of F(y) that occurs in a profile, and y ^ 
is the value of y at which it occurs. The length y is the 
normal distance from the wall and y f is the smallest value 
of y at which values from the inner and outer formula are 
equal . 


(2.36) 


(2.37) 



(2.39) 


( ^ inner 
( ^ outer 


y*y c 
y> y c 


The turbulent Prandtl number is assumed to be Pr t =0.9 
in the present study. 


2.5 Turbulence Model Implementation 

Eddy viscosity turbulence models are usually derived 
and validated for two-dimensional boundary layer flows. 
Further, the eddy viscosity coefficient determined by these 
models depends on the local flow profiles along the normal 
direction from the wall. The Baldwin-Lomax turbulence model 
performs adequately for fully attached or mildly separated 
flows over simple geometry. However, for the flows over 
more complex configurations, where the boundary layers and 
wakes may interact or flow separation may occur, the major 
difficulty encountered in applying the Baldwin-Lomax 
turbulence model in that of properly evaluating the scale 
length y ^ and in turn, of determining f° r boundary 
layer profiles. 

The turbulence length scales are determined by 1 of 
Eq. (2-30) in the inner layer, and in the outer layer. 
The eddy viscosity in the outer layer depends on the F ^ 
and the Klebanoff intermittency factor. F . l [ltr 


is a function 
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Of Ynuu and F tieb ^ a function of (y/y, ntu) • Ths y is the 
distance at which the maximum value of F(y) occurs along 
the normal direction from the wall, where F(y) is 
proportional to the moment of vorticity. For simple 
turbulence flows, a single well defined peak, exists in the 
function F(y) along a given streamwise station. When the 
flows are complex, the function F(y) may exhibit multiple 
local maxima. Selection of inappropriate length scales 
leads to inaccurate flow structure. Various methods for 
determining the appropriate length scale have been 
proposed . 

Horstman ( 1937 ) modified the Baldwin-Lomax turbulence 
model for the problem of shock-wave and turbulent boundary 
layer interaction flows. The y^ occurred outside the 
boundary layer thickness upstream and downstream of the 
shock induced separated region. The first maximum value of 
F(y) away from the wall was used to insure y^ is less than 
the boundary layer thickness. 

Degani et al.(1986, 1991) proposed modification of the 
Baldwin-Lomax turbulence model in computing the three- 
dimensional separated flow around a prolate spheroid at 
high incidence in the supersonic and subsonic flow. To 
eliminate the selection of large due to the presence of 
the vortex sheet, it was chosen at the first peak value of 
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F(y) from the wall. In case of showing a nonsmooth behavior 
in F(y) , a maximum value of F(y) was chosen at the 90% of 
the local maximum value. Another modification was that a 
cut-off distance was specified in terms of y^ from the 
previous ray. If no peak of F(y) is found in that range, 
the value of F^ and y^ t were taken as those found on the 
previous rays. 

As mentioned in introduction, the three-dimensional 
flow separation occurs in the S-duct by the pressure force 
due to the duct geometry change rather than by shear force. 
The vortical structure, which results from the flow 
separation, is stretched to the second half of the duct by 
the streamwise velocity. It causes y^ to be located outside 
the boundary layer thickness as shown in Fig. 2.1(b); 
therefore, it is not necessary to consider the whole normal 
direction from the wall to pick the correct y^. 

In this study, in order to avoid choosing an 
inappropriate length scale (y w ), the cut-off distances are 
evaluated in every crossplane. They are obtained by 
averaging the local boundary layer thicknesses within <p = 
45°. The effect of strong secondary flow due to flow 
separation can be neglected in this region. The length 
scale search is restricted to within the cut-off distance. 
If the local boundary layer thickness is less than 110% of 
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the cut-off distance, F ^ and y^ are chosen at the maximum 
peak point within that distance. Otherwise, the first peak 
value of F (y) from the wall is chosen as F m n r . 

If one employs the same method to decide a cut-off 
distance in the case of the flow with vortex generators, 
the cut-off distance is less than that of the flow without 
vortex generators. This is because the local boundary layer 
thickness within = 45° is thinner than with without vortex 
generators because the shed vortex from the vortex 
generator has a streamwise velocity deficit at the region 
of the vortex core. The cut-off distance of the flow with 
vortex generators is adjusted to be greater than the 
average boundary layer thickness which is obtained by 
averaging the local boundary layer thicknesses within <p = 


45 °. 



(a) Attached Region 


(b) Separated Region 


Fig. 2 


.1 Behavior of F(y) = y|3|[(i - exp<-yVA + ) 
the flow field 


] through 



CHAPTER 3 


NUMERICAL METHOD 


The unsteady compressible Navier-Stokes equations are 
a mixed set of parabolic-hyperbolic equations. If the 
unsteady terms are dropped from these equations, the 
resulting equations become a mixed set of elliptic- 
hyperbolic equations. These equations are more difficult to 
solve than the unsteady compressible Navier-Stokes 
equations. Most compressible Navier-Stokes equation 
solutions are obtained using the unsteady tern; the steady- 
state solutions are obtained by time marching until 
sufficient convergence is achieved. 

Both explicit and implicit schemes have been used to 
solve the compressible Navier-Stokes equations. 
MacCormack ( 1969 ) solved the compressible Navier-Stokes 
equations using an explicit scheme with a predictor- 
corrector technique. He used forward differences for all 
spatial derivatives in the predictor step while backward 
differences was used in the corrector step. Although the 
explicit schemes have an advantage that they are easy to 
implement, these schemes need long computation time because 
of the limitation on the time step due to the Courant- 
Friedrichs-Lewy (CFL) stability restriction. For this 
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reason, schemes with less-restrictive stability conditions 
have been an important subject of investigation. Allen and 
Cheng (1970) introduced a nonconsistent approximation scheme 
[/„" = (fi+i 0 -2f" +1 +fi.,“) /ax 2 ] . This approximation scheme 
becomes consistent when the steady state is reached and it 
has good stability properties when the mesh Reynolds number 
is less than 2. MacCormack(197i) modified his original 
scheme by splitting a sequence of one-dimensional 
operations. The stability condition on the revised scheme 
is less restrictive than his original scheme. Deiwert (1975) 
employed a finite volume method to solve compressible 
Navier-Stokes equations with less-restrictive stability 
condition. However, the explicit schemes are not a suitable 
method for solving high Reynolds number flows where the 
viscous regions become very thin. For these flows, a very 
fine mesh is required near the wall in order to resolve the 
boundary layer. This leads to an expensive calculation 
because of the small time step due to the stability 
restriction. 

A large and productive effort has been occurred in the 
area of implicit schemes. Polezhaev(l967) proposed an ADI 
(Alternating Direction Implicit) scheme without an 
iterative process. Briley and McDonald (1973) applied the 
generalized ADI procedure to solve the compressible Navier- 
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Stokes equation. Beam and Warming ( 1978) solved the 
compressible Navier-Stokes equation by the implicit method, 
which was the same class of ADI schemes developed by 
McDonald and Briley ( 1975 ) . MacCormack ( 198 1 ) developed an 
implicit scheme analogue of his explicit scheme. Even 
though implicit schemes are condemned for their large 
arithmetic operation counts, these schemes have been 
praised for their improved stability conditions. 

In this study, an implicit finite volume, lower-upper 
time marching code (RPLUS3D) , which was developed at NASA 
Lewis Reseach Center, is employed, and the boundary 
conditions and algebraic turbulence model are added in 
conjuction with the flow physics; the chemical reaction 
term are eliminated to save the computation time. This 
computer code employs a lower-upper (LU) factored implicit 
scheme developed by Jameson and Turkel ( 193 1) . This scheme 
is unconditionally stable in any number of space 
dimensions. Despite being implicit, the LU scheme requires 
only scalar diagonal inversions while most other implicit 
schemes require block matrix inversions. The use of scalar 
diagonal inversions offers large savings in computation 
time and temporary storage. 
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3.1 LU Scheme 

The unfactored implicit scheme with time marching for 
the vectorized Eq. (2.15) can be formulated as follows; 
viscous terms are treated explicitly to avoid complexity. 

U"* 1 = U n - a c[D K [ir^) + £ n (i v ’* 1 ) + ^(G"* 1 )] 

(3.1) 

+ aC[D ? (E") + D n (F* v ) +D { (G") ] 

where D £/ £», and D { are the spatial finite difference 

operators. The superscript n denotes the time level, i.e., 
if = U(nAt). The difficulty for solving these algebraic 
equations comes from the nonlinearity of the set of 
equation. The linearized equations with the same temporal 
accuracy can be obtain by the Taylor series expansion. 

E{ U 1 * 1 ) = E(U n ) + -2= ( C7"* 1 - U 1 * ) + 0 ( I A t I 2 ) 

dU) 

( d"* 1 - U n ) + O ( I At I 2 ) (3.2) 

n 

( u" 1 -6") + o ( I AC I 2 ) 

Let the linearized flux Jacobians of the convective flux 


Fid** 1 ) = F{U n ) + 21 

[ dU 


= G(l 7°) + 

l dU 


vectors be 
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A = 


a| 

dU 




(3.3) 


and the correction be A U = U D+1 - U “. The Jacobian 
matrices A, B, and C are given in Appendix A. 

This unfactored implicit scheme is first order 
accurate in time. Therefore, the second higher order terms 
can be neglected without loss of time accuracy of the 
linearized governing equations. 


[J + 4t(D E A + D n B * D( C ) ]aZ7 
= -a t[D i (S-E v ) * (F-F v ) +D : (G-G V )} 


(3.4) 


The linearized Eq. (3.4) with the unfactored implicit 
scheme has large block banded matrices, which require large 
storage and computation time for inversion. The Eq. (3.4) 
can be factorized by replacing the operator with a product 
of three one-dimensional operators. This is sane as the ADI 
scheme, which also requires inversions of block tridigonal 
or block pentadiagonal matrices. If one solves the block 
tridiagonal system by Gaussian elimination without 
pivoting, the operation count for the block Thomas 
algorithm is 0(NM i ) where M is the block size and U is the 
number of unknowns. Clearly it is desirable to avoid 
solving a block tridiagonal system. For many standard 
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algorithms, one can not be confident that the Thomas 
algorithm is numerically stable if the diagonal dominance 
is lost by increasing time steps. 

Jameson and Turkel(1981) proposed the idea of a lower- 
upper factored implicit scheme that is unconditionally 
stable in any number of space dimensions and also yields a 
steady-state solution that is independent of At. The LU 
implicit scheme needs only two factors even for three- 
dimensional problems because of the unique manner of 
splitting. As a result, this scheme is more stable and 
robust than ADI schemes. Let 

A = A * + A 

B = B' + B~ (3.5) 

C = d* + C' 

The split flux Jacobians, A + , B + ,C + , A', B' and C* are 
constructed such that the eigenvalues of "+" matrices are 
nonnegative and those of matrices are nonpositive. Of 
the many ways of splitting, Jameson and Yoon's(1987) method 
is employed as follows: 
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A" = 0.5 (A + y x I) 

A' = 0.5 (A - y x I) 

B* = 0.5 ( B + y§I) 

(3.6) 

S' = 0.5 ( B - y §1 ) 

C’ = 0.5 ( C + y d J) 

C' = 0.5 ( c - YeJ) 

where 7 ^, 7 * and 7 C are greater than the spectral radii of 
the flux Jacobians associated with them: 

y & = max ( | X A | ) 

Yj -max ( |X,|> (3 . 7) 

Yc = ^ax ( | *e| ) 

Here, X^, X s and X^. represent eigenvalues of Jacobian 
matrices A, B and C. The spectral radii and eigenvalues of 
Jacobian matrices are obtained in section 3.3. 

Substituting Eq. (3.5) into Eq. (3.4) and performing 
the first order upwind difference according to the sign of 
the eigenvalues, The linearized implicit scheme can be 
obtained : 

[ J + A C ( D{A* + D'B* ♦ D{ C * * D{A~ + D*B' + D{C~ ) ] a U 

= -ot:[D l (E- E v ) + D n { F - F v ) +D Z (G-G V )} (3.8) 
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where £> { ‘ , , and d: denotes backward-difference operators 
while D { + , , DS are forward-difference operators. Eq. 
(3.8) can be expanded in discretized form as follows: 


♦ At (A,v.*A 0 iJJC - A-i^ACJ..^ + 

A* ljjtA&i.tjj 

^AU, v ^B^AE7. 

/.* A - lj.k ijjc + 

JW&W 


i.k ~ lj.k AJ,. + 

\j.k A J, . t 

- 

. AU lVJt ) = A t J?HS 



(3.9) 


where 


RHS=D i {E-E v ) + D n (F - F,) * D f ( G - G J (3.10) 

This discretized equation can be written as 

[ { X + At (A* -A" + B * - B~ + C* - C* ) } 

f At ( D £ 'A* + B' + n: C' - A* - B* - C* ) (3.11) 

+ At ( £* A' + D„* B' + D r * C‘ + A* + B' + C* ) ] A£7 

= AtRHS 


Eq. (3.11) can be factorized according to the sign of the 
Jacobian matrices 

{X + At ( D t " A * + D,' B * + £>.' C* - A* - B* - C* ) } ( JC ) ■ 

{JC + At ( A* + D,*B‘ + Df C* + A* + B* + C* ) }A£7 


where 


JT = X + At ( A* - A* + B* - B~ + C* - C' ) (3.13) 


Notes that matrix is diagonal. This can be easily 
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verified by substituting Eq. (3.6) into Eq. (3.13). K ] is 
also diagonal and can be moved to the right hand side. 

{ I + At ( D" A* + + D r 'C* - A' - S' - C* ) } 

{l + At ( D t * A" + D*J§ _ + D r * C* + A' + B* + C* ) }AC7 

= At {J v At (y A * y 6 + 7 £ ) X (3>14) 

The operator on the left hand side of Eq. (3.14) 
represents Lower and Upper operator of this scheme. These 
two operator represent forward and backward substitutions. 
It is interesting to note that if there is no source term 
in the governing equation, the numerical method completely 
eliminates the need for block matrix inversion. In fact, 
the two operators in Eq. (3.14) require only scalar 
inversions. Although the LU scheme is an implicit scheme, 
the numerical operation counts are not much differrent from 
those of explicit methods. 

The discretized equation in the finite volume method 
is derived by approximating the integral form of the 
equation to be solved. The computational region is divided 
into elementary quadrilateral volumes within which the 
integration is carried out, and the integral equation is 
evaluated at each subdomain. This method can easily handle 
the complicated geometry without considering the equation 
written in curvilinear coordinates. It makes it possible to 
avoid problems with metric singularities that are usually 
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associated with finite difference methods. If this method 
is applied on the uniform rectangular grid, the discretized 
equation will be equal to the discretized equation using 
the central finite difference method. It has second order 
accuracy in space, but for the non-uniform grid the 
convergence rate in space is less than second order. 


3.2 Artificial Dissipation 

The finite volume formulation reduces to a central 
difference approximation on a uniform grid. It allows 
undamped oscillations with alternative signs at odd and 
even mesh points. Wiggles appear in the neighborhood of 
severe pressure gradient regions or stagnation points. 
These spurious oscillations can not be smoothed out totally 
by the viscous and dissipation terms. In order to suppress 
these numerical oscillations, the artificial dissipation 
terms are added into the LU scheme. 

In this study, Jameson' s ( 198 1) adaptive artificial 
dissipation scheme is employed. The dissipation terms 
consist of blended second and fourth order differences. The 
fourth order difference terms provide background 
dissipation throughout the flow field to prevent odd-even 
decoupling which occurs from the linearized Euler equation 
terms. The second order dissipation terms are used to 
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stabilize the flow calculation near the regions of the 
strong pressure gradients. These terms are explicitly added 
to the RHS term as an additional residual. The added 
dissipation term are as follows: 


* < «"> Cj.,. A * *<§> i. ». * > 


(3.15) 


where 


e' 21 ! 




v i«a, j , k ' v i , j , 


v i- i.j.k) 


(3.16) 


1» J, * 


^ Pi, j, Jc * -Pj-l, j, Jc 
Pi-1, j.k + ^Pi.j.k* Pi-1, j, Jc 


(3.17) 




= min ( ( J"? ) k ) 


(3.18) 


€ 



J. * 


max ( 0 , k. - e 121 ! . ) 

4 i--.j.k 


(3.19) 


and k 4 are scalar constants. In this study, k 2 and * 4 are 
1 and 1/32, respectively, and the magnitudes of artificial 
dissipation coefficients are much less than the eddy 
viscosity in the boundary layer. The term -y iik is a spectral 
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radius scaling factor and is defined as 

h.j.k = + Yj + Ye 

= |5 jr ||u| + |5yllv| + |5 I lk| + c(U + 5j + 5i>^ 

x (3.20) 

+ hxl M + h/l M + hzl I W \ + C ( nx + T)y + ) 2 

+ |CJ|U| - |C y | I V| + |C,| \w\* c ( Cj* Ci ) * 

which is the sum of the spectral radii of A, B and C. 

The first terms in the parentheses of Eq. (3.15) are 
the second order dissipation. It has an extra pressure 
gradient coefficient which is constructed by taking the 
second difference of the pressure. Its value increases in 
the neighborhood of the strong pressure gradient region, so 
the non-physical overshoot or undershoot are eliminated by 
the second order term. The second terms in the parentheses 
of Eq. (3.15) are the fourth order dissipation. The 
coefficient £ (4) switches off when the second order nonlinear 
coefficient is larger than the constant of the fourth order 
coefficient. 


3.3 Eigenvalues of Jacobian Matrices 

The eigenvalues of Jacobian matrices are required to 
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analyze the stability of a numerical scheme. The Jacobian 
matrices of non-conservative equation (A, B and C) are much 
simpler than the Jacobian matrices of conservative equation 
(A, B and C) . Warming et al.(1975) showed that an uniformly 
bounded similarity transformation between the Jacobian 
matrices of non-conservative equation and conservative 
equation existed in the invicid gas dynamic equations. The 
Jacobians of the generalized trnsformed convective flux 
vectors can be expressed by the Jacobians of conservative 
equation . 


C 3 Z* _ r 3 Ej r 3 F* r 3G _ j , r Dj_P /— 

A = Jo - 5 x ~dU '*Ju * ,Z ~3U ~ 




A 3G r BE t dF r 3G r ■» r r * 

c - Tu - W * t'TTu ' V 


Using the similarity transformation, Eq. (3.21) is 
changed to the simple form. 

A = M ( * Zy3+ l z C) M' 1 

B = M {'r\ x A + r\j3 + r\ z C) M' 1 (3.22) 

C = AT ( + CyS + C Z C) AT 1 


It is not hard to find the eigenvalues of the Jacobian 
matrices in the non-conservative equation. The eigenvalues 
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°f A are as follows: 

^•1,2,3 = u 

1 4 = U + C (3.23) 

1 5 = u - c 

The eigenvalues of B and C are similar to those of A, 
only u in the eigenvalue of A has to replace to v and w, 
respectively. 

The eigenvalues of A, B and C are easily obtained 
using the Eq. (3.23). The eigenvalue of A are as follows: 

* 1 . 2.3 = + *-y v + Zz w 

X 4 = 5 X U + l y V+ \ z w+ >* (3 *24) 

*5 = *>x U + ty V + Zz W ~ C ( + + ) 2 

The eigenvalues of B and C are similar to those of A, 
only £ in the eigenvalues of A has to replace to rj and 
respectively . 

In order to obtain the bound of the spectral radii in 
the flux Jacobians, the biggest eigenvalue is tested. 
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Yi = max ( X i ) i=ljS 

* k j i u\ * ig i v| + 15,| i c ( e x - e y + ) * 

£ AAA ( | u i + | v| + | w\ +/3c) 

£ AAA ( 2 \/u 2 + v 2 + w 2 +2 c ) 


where AAA = max(£ x , £ y , . 

Using the same method, the bound of the spectral radii 
of B and C can be obtained: 

y § = BB3 ( 2 i/u 2 + v 2 + w 2 +2 C ) (3.26) 

= CCC( 2 Vu 2 + V 2 + w 2 + 2 c) (3.27) 

where BBB = max(i 7 x , 7 j y , 77 J and CCC = max(f x , f y , 
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DIFFUSING S-DUCT WITHOUT VORTEX GENERATORS 


4.1 Geometry and Grid 

The geometry of the diffusing S-duct examined in this 
study is shown in Fig. 4.1. The duct centerline is defined 
by two circular arcs with identical radii of curvature, 
which are 5 times the inlet duct diameter, and subtended 
angle 6 m „/ 2 = 30°. Both arcs lie within the xy-plane as 
shown in Fig. 4.1. The coordinates (x cl , y cl , z cl ) of the duct 
centerline are given by Eq. (4.1): 

For 0 < 6 < 0 mx /2 

x c , = R sin 0 
y cl = R cos 0 - R 
Zd = o 

For e^l 2 < 9 < (4.1) 

x ci = 2 R sin ( ) - R sin ( 0^ - 0 ) 

0 

y c i = 2Ecos ( -S5 ) - R - Rc os (0^-0) 

2 ci = 0 

The cross-sectional shape of the duct perpendicular to 
the centerline is circular. The diameter of the cross 
section varies with the arc angle 6 and is given by Eq. 

(4.2) . 
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D 

£> 


= 1 * 3 ( f- 1) < r 


2 { Tr 1) ( <C 


(4.2) 


where D ; and D e are the diameter at the S-duct inlet and 
exit, respectively. The area ratio of the duct exit to 
inlet is 1.51. The offset of the duct resulting from the 
centerline curvature is 1.34D ; . The length of the duct 
measured along the centerline is 5.24Dj. A straight pipe, 
which is 4.6D, long, is installed upstream of the S-duct to 
provide the desired boundary layer thickness at the inlet 
of the S-duct. In order to minimize any downstream effect, 
a 9D e straight section of pipe is attached at the exit of 
the S-duct. The average inlet Mach number is 0.6 and the 
Reynolds number based on the duct diameter is 1.76xl0 6 . 

In the present study, an 0-grid is adopted because it 
conforms well to the boundaries of the circular duct. The 
O-grid consists of 47 radial points, 42 circumferential 
points in the half duct, and 70 streamwise nodal points. A 
finer grid is used in the region of flow separation. 
Exponential stretching is used to obtain a fine mesh near 
the wall. The upstream and downstream lengths of straight 
ducts are also extended using the exponential stretching. 
The mesh size adjacent to the duct surface is almost 
1.25X10" 4 times the duct inlet diameter. The two grid points 
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nearest the wall are at value of y + of about 2.6 and 5.7 at 
the reference station (S/Dj = -1.5). 

The computed results do not depend on the initial 
velocity conditions, i.e., the initial velocity profile 
with or without adjusting the axial velocity by the one- 
seventh power velocity distribution law near the wall. The 
mass flow changes between the inlet and exit was within 1 
percent for all calculations. The residuals for these 
numerical solutions were reduced approximately three orders 
of magnitude. Solutions were obtained on the Cray-YMP. The 
number of iterations required to obtain the converged 
solutions was approximately 40,000. The computational speed 
was approximately 960 iterations per CPU hour. 


4.2 Results and Discussion 

When discussing numerical and experimental results, 
streamwise position will refer to the distance to cross 
stream-planes measured from the inlet of S-duct along the 
duct centerline and normalized by the duct inlet diameter. 
Position within cross stream-planes is specified by the 
polar angle <p, measured from the vertical in a positive 
clockwise direction as shown in Fig. 4.1, and the radial 
distance from the centerline of the duct. 
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Fig. 4.2 shows the surface static-pressure 
distributions at <p = 10° , 90° and 170° which are compared 
with two experimental data. Note that the definition of the 
surface static-pressure coefficient in the two experiments 
is different. Vakili et al.(1986) measured the reference 
flow parameters at S/D ; = -1.5, upstream of the S-duct for 
normalizing downstream flow data. The reference variables 
were evaluated at the center of the duct. 

Wellborn et al.(1992) measured the reference flow 
parameters at S/D ; = -0.5, upstream of the S-duct. The 
reference dynamic pressure was evaluated by subtracting the 
wall static-pressure from the total-pressure measured at 
the center of the duct. They used a similar duct but larger 
than that used by Vakili et al.(1986); therefore, the 
Reynolds number of Wellborn et al. (1992) experiment is 47% 
higher than that of Vakili et al.(1986) experiment. 
However, In this study, calculations were made using the 
same Reynolds number as the Vakili et al.(1986) experiment. 

The computed surface static-pressure distributions are 
in good agreement with the experimental data except in the 
separation region. In the separation region, the predicted 
values of surface static-pressure are higher than the two 
experimental results. Both experimental data show constant 
values of static-pressure at o - 90° and 170° in the region 
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2 < S/D; < 3; the computational result shows a similar 
result in the region 3 < S/D, < 4. 

The experimental flow separation region shown in Fig. 
4.2 was determined by surface oil flow visualization. The 
computed flow separation region is determined by examining 
the streamwise velocity in the vertical plane of symmetry. 
The predicted separation length is 1.94, which is a little 
shorter than the experimental value of 2.1. The predicted 
separation (2.44 < S/D; < 4.40) occurs farther downstream 
than was observed experimentally (2.02 < S/D; < 4.13). This 
indicates that the applied turbulence model, even as 
modified, cannot correctly account for the three- 
dimensional separation flow with very strong secondary 
flow. The experimental and numerical results show that the 
flow fields in a diffusing S-duct have strong secondary 
velocities with flow separation, and the counter-rotating 
vortices resulting from the flow separation are stretched 
into the second half bend of the duct by the streamwise 
velocity. These complex flow fields result in the moment of 
vorticity (F(y) ) having several peak values along the normal 
direction from the wall. Although the first peak value from 
the wall is chosen as the length scale (y,^) in order to 
avoid choosing an inappropriate length scale, this chosen 
length scale in the flow separation region cannot be 
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considered as a perfectly correct length scale. 

The reverse flow in a diffusing S-duct is associated 
with the adverse pressure gradient due to the increase of 
duct area and the secondary flow due to the duct curvature. 
Fig. 4.3 shows the velocity profiles in the vertical plane 
of symmetry. The reverse flow occurs away from the wall; an 
enlarged view is shown in Fig 4.3(b) to display this 
feature more clearly. These different characteristics of 
flow separation can occur due to the turbulence model. If 
the function F(y) has a peak value close to the wall, the 
eddy viscosity along the normal direction from the wall 
approaches quickly to zero by the Klebanoff intermittency 
factor except the region of the near wall. This incorrect 
viscosity profile cannnot adequately account for the 
reverse flow associated with the adverse pressure gradient 
and the strong secondary flow. 

Fig. 4.4 shows the surface static-pressure 
distribution along the circumferential direction at the 
three different streamwise locations S/D; = 0.96, 2.97 and 
4.01. The computational results at S/D; = 0.96 and 4.01 
agree quite well with the experimental data. S/D ; = 0.96 and 
4.01 are located upstream and downstream of the flow 
separation region, respectively. The computed values of 
surface static-pressure at S/D; = 2.97, which is located 
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within the flow separation region, are higher than the 
experimental data. This overprediction of surface static- 
pressure seems to result from the inadequate turbulence 
model as previously mentioned. 

Fig. 4.5 shows the static-pressure contours at the 
various streamwise locations. The computed results are 
compared with Vakili et al.s(1987) experimental data. Since 
the flow is symmetric with respect to a vertical plane 
passing through the centerline, only half of the plots are 
shown in these figures. The calculated static-pressure 
contours show similar trends as the experimental results, 
but the computed static-pressure levels are higher than the 
experimental values. The static-pressure coefficient are 
evaluated as (p iocaJ - p rcf ) /q ref , and the reference values are 
measured at the center of duct in the reference plane (S/D ; 
= -1.5). Comparing two experimental results of the surface 
static-pressure coefficient of Fig. 4.2(a) and the static- 
pressure coefficient contours of Fig. 4.5, the static- 
pressure coefficient near the wall in Fig. 4.5 is much 
lower than that shown in Fig 4.2(a). However, the computed 
static-pressure coefficient near the wall in Fig. 4.5 are 
very close to the experimental surface static-pressure 
coefficient, and also Fig. 4.2(a) shows that the surface 
static-pressure coefficients, even if at the reference 
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plane, are much different along the circumferential 
direction. This probably results from deficiencies in the 
experiments, primarily coarse data acquisition locations 
and uncertainties in the static-pressure measurements using 
pitot tubes. 

Figs. 4.5(a) and 4.5(b) show the increase of the 
static-pressure toward the outer wall in the first half 
bend. This result is anticipated by the inviscid theory. In 
the second bend, the static-pressure increases from the 
upper wall to the lower wall as shown in Figs. 4.5(d) and 
4.5(e) due to the adverse curvature direction. The static 
pressure along the duct increases due to the increase of 
duct area. The static pressure core shown in Fig. 4.5(e) 
results from the streamwise velocity deficit at the region 
of the two counter-rotating vortices. This means that 
nonuniform flow at the exit occurs from the flow 
separation . 

Total-pressure contours compared with the experimental 
data obtained by Vakili et al.(1987) are shown in Fig. 4.6. 
Fig. 4.7 shows the total-pressure contours compared with 
the experimental data obtained by Wellborn et al.(1992). 
Note the different definition of the total-pressure 
coefficient in the two experimental data. The agreement 
between the computational and experimental results is quite 
good except downstream of the flow separation. The 
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disagreement at the downstream of the flow separation 
caused by the different flow separation region. 

A qualitative picture of the secondary flow pattern 
in a curved duct is that an inviscid core fluid moves 
toward the outer wall of the duct, and a low speed boundary 
layer fluid migrates circumferentially from the outer wall 
to the inner wall in the first half of the S-duct. This 
phenomenon results in low energy flow accumulating near the 
inner wall of the first half bend. This is shown in the 
total-pressure contours of Fig. 4.6(d). The adverse 
pressure gradient is induced on the second half bend of the 
duct due to increase of the duct area. The pressure 
gradient causes a thick boundary layer and deflection of 
the streamwise flow direction. 

The above mentioned secondary flow pattern contributes 
to the formation of a pair of counter-rotating vortices by 
the three-dimensional flow separation. Tobak and 
Peake (1982) showed the topographical structure of three- 
dimensional flow separation. The counter-rotating vortices 
formed by the vortex lift-off stretch to the exit of the S- 
duct by the streamwise velocity, and move away from the 
wall to the center of the duct. In the region between two 
counter-rotating vortices, the secondary velocities induced 
by these vortices push the low energy flow toward the 
center of the duct. The high energy flow between the 
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vortices and the duct wall is pushed toward the boundary 
layer. This mechanism makes the convex shape of the 
inviscid core flow region as shown in Fig. 4.6(f). 

The shape of the total-pressure contours in the cross 
plane depends on the strength of the counter-rotating 
vortices and the core location of the vortices in that 
plane. They are related to the original location of the 
counter-rotating vortices in the duct. The computed three- 
dimensional flow separation region occurred further 
downstream than was observed in the two experiments. This 
causes the discrepancy between the computational and 
experimental total pressure contours at S/D; = 5.24 and 
5.73 . 

Comparing Figs. 4.7 and 4.8, we see that axial Mach 
number contours are very similar to the total-pressure 
contours at the same axial location. The computational 
total-pressure contours at S/D ; = 5.24 and 5.73 indicate 
that the computed streamwise velocity deficit (U,* - u) at 
the region of the counter-rotating vortices is bigger than 
was observed experimentally. This large streamwise velocity 
deficit makes the inviscid core flow region larger in order 
to satisfy the constant mass flux along the streamwise 
direction . 

Fig. 4.9 shows the secondary velocity profiles at the 
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five stations along the duct. They are compared with the 
experimental results obtained by Vakili et al.(1987). Fig. 
4.10 is the secondary velocity profiles compared with the 
experimental results obtained by Wellborn et al.(l992) at 
S/Dj = 5.73. The development of secondary flow in the curved 
duct is clearly shown in these figures. The computational 
results are in good agreement with experimental data except 
downstream of the flow separation region. The secondary 
velocity profiles in the first bend clearly depict the 
qualitative picture of the secondary flow pattern in the 
curved duct as mentioned in the discussion concerning the 
total-pressure contours. 

Fig. 4.9(c) shows the accumulation of low energy flow 
at the lower wall, which is consistent with the observation 
of the total-pressure contours. Downstream of the flow 
separation region, Figs. 4.9 and 4.10 show that a pair of 
counter-rotating vortices move away from the wall and 
toward the center of the duct. The computational results 
show that the secondary velocity is overestimated 
downstream of flow separation. This results from the small 
eddy viscosity effect in the flow separation region by the 
implemented turbulence model, i.e., F ^ and are chosen 
at the point of the first peak value from the wall in that 
region. 
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The variations of boundary layer thickness at <p = 10°, 
90° and 170° along the duct are shown in Fig. 4.11. The 
boundary layer thickness is defined as the normal distance 
from the wall where the total-pressure coefficient is 1.0. 
The predicted results are compared with the experimental 
results obtained by Vakili et al.(19S7). The computational 
results and experimental data are in reasonable agreement. 
The rapid boundary layer growth at <p = 17 0° is caused by 
the flow separation. In the transition region (S/D; = 0.0) 
from the straight duct into the first bend, the computed 
results show that the boundary layer thickness at <p = 170° 
is less than that at ^ = 10°. The streamwise velocity near 
the lower wall in the transition region is faster than that 
near the upper wall due to the effect of the curved 
geometry. It was well depicted in the static-pressure 
contours as shown in Fig. 4.5. The experimental data do not 
clearly show the effect of this flow mechanism. As shown in 
the secondary flow pattern of Fig. 4.9, the high energy 
flow migrates toward the outer wall in the first bend, 
therefore the boundary layer thickness at = 10° along the 
duct is less than that at ip = 90° and 170°. 

Downstream of the flow separation, the computational 
result shows that the boundary layer thickness at <p = 90° 
is less than that at = 10°. The reason is that the strong 
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secondary velocities induced by the counter-rotating 
vortices push the high energy flow toward the wall. The 
stronger secondary velocities, as compared with 
experiments, are associated with the fact the computed 
total-pressure boundary layer is thinner than the 
experimentally measured one. 

Fig. 4.12 shows the velocity profile along the normal 
direction from the wall at the four different streamwise 
locations. The first grid points in the computation are 
located inside the visccus sublayer (y + < 5) . At the first 
grid points, the friction velocity is calculated to 
normalize the velocity profile. The viscous sublayer 
region, the log linear region, and the wake region are 
shown in this figure. In Fig. 4.12(b), the velocity profile 
at <p = 170° is not shown because the definition of the 
friction velocity is not applicable in the flow separation 
region. At the exit of S-duct (S/D ; = 5.2), the flow is 
reattached but a pair of counter-rotating vortices are 
present as shown in Fig. 4.9(e). These cause the boundary 
layer profile to deviate from the law of the wall at <p = 
170°. The velocity profiles at <p = 170° show the large 
streamwise velocity deficit (U a - u) . Fig. 4.12(d) shows a 
comparison with the velocity profile measured by Wellborn 
et al.(1992) at S/D- = 5.73. The agreement in the wake 



67 


region is poor because the strength of the counter-rotating 
vortices was overestimated as previously mentioned. 

The skin friction values along the streamwise 
direction are plotted for <p = 10°, 90° and 170° in Fig. 
4.13. Note that there is no experimental data for the skin 
friction values. The trends of the computed results are 
similar to the trends of Bansod and Bradshaw's(1973) 
experimental data for low speed flow in a nondiffusing S- 
duct . The skin friction decreases along the duct due to the 
increase of the duct area. 

Fig. 4.14 shows the streamlines in the symmetry plane 
along the duct. The experimental result was obtained by 
placing a thin metal plate in the symmetry plane of the S- 
duct. Even though there is no cross flow in this symmetry 
plane, the presence of thin plate in the symmetry plane 
introduces shear layer development and blockage. However, 
the comparison with Wellborn et al. '3(1992) experimental 
result agrees well qualitatively. 

In the current computations, numerical results 
demonstrate the capability of a modified algebraic 
turbulence model in the flow fields of the three- 
dimensional flow separation with a strong secondary flow. 
The computed results agree quite well with the experimental 
results except in the flow separation region. Even though 
there are deviations between experimental and numerical 
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results in the flow separation region, the computed results 
depict well the flow structure in the diffusing S-duct. 
However, further studies to obtain the correct length scale 
in the flow separation region are required. 
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Fig. 4.2(a) Axial surface-static pressure coefficient. 

Cp, = <Plocal - P^) 

( Ma = 0.6, Re d = 1.7 6x10* ) 

Exp. Vakili et al. (1986) 
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Fig. 4.2(b) Axial surface-static pressure coefficient. 

^p2 (Plocai Pwafl) / ^P«^l ” P wail ) 

( Ma = 0.6, Re d = 1.76x10* ) 


Exp. Wellborn et al.(i992) 

( Ma = 0.6, Re d = 2.6x10* ) 



(a) at Vertical Plane 



Fig. 4.3 Streamwise velocity profiles in the vertical 
plane of symmetry on the S-duct. 

( Ma = 0.6, Re d = 1.76x10* ) 
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.8 Axial Mach number contours without vortex 
generators. At S/D ; = 5.73 
( Ma = 0.6, Re d = 1.7 6X10 4 ) 

Exp. Wellborn et al.(l992) 

( Ma = 0.6, Re d = 2.6X10 4 ) 
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Fig. 4.10 Secondary velocity profiles without vortex 
generators. At S/D, = 5.73 
( Ma = 0.6, Re d = 1.76x10 s ) 

Exp. Wellborn et al. (1992) 

( Ma = 0.6, Re d = 2.6X10 4 ) 
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Fig. 4.11 Boundary layer thickness. 

( Ma = 0.6, Re d = 1.76x10* ) 

Exp. Vakili et al. (1987) 









CHAPTER 5 


STRAIGHT DUCT WITH VORTEX GENERATORS 


The spiral longitudinal vortex interactions with 
turbulent boundary layer are numerically investigated in a 
cylindrical duct. The helical motion of the injected vortex 
is compared with the prediction by imagine vortex system 
and the prediction by Wendt et al. / s(1992) vortex 
interaction model. Two prediction models are derived in the 
Appendix B and C. In a second model, the constants which 
were derived from the experimental result of the external 
flow are employed. Although it is not sufficient to apply 
the same constants to predict the helical motion of the 
injected vortices in the internal flow, a reasonable 
prediction can be obtained in a short region just 
downstream of the vortex generators. 

Kunik(1986) conducted a numerical study about the 
behavior of the injected vortex using the PNS equations on 


the straight duct. The 

flow was 

incompressible 

and 

the 

Reynolds 

number 

based 

on pipe 

diameter 

was 2000. 

The 

injected 

vortex 

was 

set up 

at the 

inlet 

of 

the 


computational region because the PNS equations could be 
solved by forward marching in space. Note that the PNS 
equations cannot consider the streamwise velocity deficit 
at the vortex core sufficiently because of neglecting the 
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streamwise diffusion term of the FNS equations. 

In the present study, the injected vortices are set up 
within the computational region, and three-dimensional FNS 
equations are solved by the previous described numerical 
technique. Fig. 5.1 shows the computational grid for a 
cylindrical duct with L/D = 20.0. The polar grid topology 
consists of 47 radial points, 73 circumferential points and 
60 streamwise nodal points. Exponential stretching is used 
to obtain a fine mesh near the wall. In order to obtain 
high quality velocity profile, the wall shear stresses are 
measured within the viscous sublayer. The first grid point 
nearest the wall has a y + value of less than 3, which is 
about 1.6x10"* times the duct diameter. The location of 
vortex generator is at X/D = 2.1. The entrance Mach number 
is 0.6 and the Reynolds number based on the inlet diameter 
is l.OxlO 6 . 

The number of iterations required to obtain a 
converged solution was approximately 25,000. Solutions were 
obtained on the Cray-YMP. The computational speed for the 
full duct was approximately 540 iterations per CPU hour. 
The residuals for these solutions were reduced by almost 
three orders of magnitude. The mass flow changes between 
the inlet and exit were within 1 percent. 
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5.1 Vortex Generator Model 

The shed vortex from the vortex generator is modeled 
by providing the two-dimensional secondary flow structure 
in a crossplane. The secondary velocity structure is 
formulated as a viscous trailing vortex on the assumption 
of steady, incompressible, laminar and axisymmetric flow. 
The secondary velocity structure obtained with the above 
assumptions can be applied to the compressible and 
turbulent flow, because only one crossplane of the 
computational domain employs this vortical structure to 
simulate the shed vortex downstream of vortex generator. 

The Navier-Stokes equations in cylindrical coordinates 
based on the origin of the trailing vortex in the infinite 
space are as follows: 


,. , 3u r du, 

radial worn. eg. u r + u r - 

r dr ox 


+ v [ V 2 u r 



1_ Bp 
p dr 

(5.1) 


rotational mom. eg. 


u. 


0U 9 

dr 


+ u. 


3uq 

dx 


r 


= V [ V 2 U 9 



] 


(5.2) 



axial mom. eg. 


du 

X 

dr 
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du x 

dx 


1 dp 
p dx 


+ v V 2 u x 

(5.3) 


1 d(ru.) du x 

continuty eg. = — — + = 0 

r dr dx 


(5.4) 


where 


V* . JL . _J E. 

dr 2 tdr dx 2 


(5.5) 


These equations are linearized and solved by making 
the following assumptions: 

1) The streamwise velocity deficit u d = U m - u x and the 
rotational velocity u 9 are small compared to the 
free-stream velocity U a . 

2) The radial velocity u r is very small compared to U„. 

3) The Reynolds number of the main flow, U^x/v, is 
large . 

These assumptions reduce the above momentum and continuity 
Eqs . (5.1) - (5.4) to 


radial mcm. eg. 


U 8 Z _ _1 dp 
r p dr 


(5.6) 


rotational mom. eg. 


U m 


dUg 

dx 


v ( 


d 2 ^ 

dr 2 


1 3Ug Ug 

r dr r 2 

(5.7) 
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axial mom. eg. 



8uj ^ i 9u d . 
dr 2 r dr 


(5.8) 


1 d(ru r ) du d 

continuty eg. .. ■ - -r-= = 0 

r dr dx 


(5.9) 


The boundary conditions to be satisfied by these 
linearized equations are: 

x > 0 , u B - 0 and u d - 0 as r - « 

x - <» , u 9 — 0 and u d — 0 for all r 

x = 0 , u, = T/2nr , u d = 0 except at the singular 

point r = 0. 

By the nature of the approximations, the vortex is examined 
at some distance downstream of its origin. Hence it is 
sufficient to assume that the vortex is suddenly generated 
at x = 0 as a free vortex of circulation r. Far downstream, 
the vortex finally decays until all the perturbation 
velocities u r , u ( and u d are once again zero. Under these 
boundary conditions, the solutions of the reduced Eqs. 
(5.6) - (5.9) are as follows: 


2nr 


[ 1 - exp ( 


•^-r : 

4VX 


) ] 


(5.10) 


u 


r 


Ar 

2x 2 


exp ( - 


U.r‘ 


) 


4vx 


(5.11) 
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u d = - exp ( - 


C/.r 2 

4vx 


(5.12) 


The integration constant A can be found from equating the 
change of momentum of the flow in the entire wake flow to 
the drag on the vortex generator; 

A = — — (5.13) 
4upv 


where D a is the profile drag of the vortex generator. 

Rotational velocity Eq. (5.10) and radial velocity Eq. 
(5.11) can be used to set up the vortical flow in the cross 
plane. However, comparing the magnitude of these velocities 
with the assumptions that a small section of NACA 0012 wing 
is used as a vortex generator with a proper angle of attack 
and the length x from the origin of vortex to the cross 
plane for vortical structure is 0(1), the radial velocity 
is small compared with the rotational velocity. 



U r 2 

[ 1 - exp ( - — = — ) ] 
4vx 

Ar , U„r 2 

exp ( - — ) 




£££2: 1 1 . _L , ) . J : { Ezll 

D 0 21 4vx 3 ! 4VX 
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Hr 2 

4vx 



t/.r 2 

4vx 


+ . • • ] * O(10 2 ) 

(5.14) 


Only the rotational velocity is used to make the vortical 
structure in a crossplane. 


u e 


r 

2 m r 


[ 1 - exp ( - 


U m r 2 

4vx 


) ] 


(5.10) 


If we apply the vortical structure of Eq. (5.10), 
which is formulated as one fully rolled up trailing vortex, 
to the circular duct directly, normal velocity component 
exists on the duct wall. In order to consider the shed 
vortex created from the vortex generator mounted within the 
circular duct, we can employ the image vortex because of a 
very small vortex core just downstream of vortex generator. 
The image vortex of equal strength as the inviscid flow is 
located outside the duct using Milne-Thomson ' s circle 
theorem (1968) . The superposed vortical flow within the duct 
has no normal velocity component at the wall. The 
tangential velocity component approaches zero at the wall 
by reducing the magnitude of the superposed vortical flow 
inside the boundary layer by the one-seventh power velocity 
distribution law. These adjusted vortical velocities are 
introduced at every point in the crossplane. 

In order to consider the streamwise velocity deficit 
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(u^ = U„ - u x ) in that crossplane, the flow passing the 
vortex generators is assumed as the steady-state steady- 
flow process. The temperature and stagnation enthalpy at 
the crossplane are calculated by averaging these values in 
the upstream-plane and downstream-plane of the crossplane 
for vortical structure. Even though the streamlines between 
two crossplanes are not the same as the streamwise 
direction due to the vortical flow, this approximation is 
sufficient if the vortical flow is small compared with the 
axial flow or the distance between two planes is small 
compared with the duct diameter. The stagnation enthalpy 
obtained by this approximation are uniform at every local 
grid points of the crossplane because of the streamwise 
velocity deficit in the downstream-plane for vortical 
structure. From this stagnation enthalpy at the local grid 
point, we can obtain the deficit of streamwise velocities 
in the crossplane with the calculated vortical velocities 
and temperatures. 


5.2 Results and Discussion 

In order to examine the usefulness of the vortex 
generator model and to investigate the effect of the 
different type vortex generators in a straight duct, four 
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different cases are tested: (1) single embedded vortex, (2) 
counter-rotating vortices of the same strength that rotate 
toward each other, (3) counter-rotating vortices of the 
same strength that rotate away from each other, and (4) co- 
rotating vortices, one vortex having double the strength of 
the other. 


number 

T/DU. 


X = r/R 

(1) 

0.062 

0° 

0.831 

(2) 

-0.062, +0.062 

i 

-18°, 13° 

0.831 

i 

(3) 

1 

+0.062,-0.062 

! 

1 -27°, 27° 

0.831 

(4) 

-0.062,-0.031 

| 

| -27°, 27° 

0.831 


Table 5.1 The strength and location of the embedded 
vortices (T is positive when the vortex rotates counter- 
clockwise, and tp is the circumferential angle from the 
vertical plane on the lower wall) 


The boundary layer thickness at the axial location of 
the vortex generator (x/D = 2.1) is 0.06 times the duct 
radius. The vortex generator is at a height of 0.16 times 
the duct radius. Therefore, the vortex generator tip is 
located well outside of the boundary layer. 


Figs. 5.2 - 5.4 show the computational results when a 
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single vortex is embedded in a crossplane within the duct. 
The total-pressure contours and secondary velocity profiles 
at the several different streamwise locations are shown in 
Fig. 5.2 and Fig. 5.3, respectively. Fig. 5.2(a) and Fig. 
5.3(a) are the total-pressure contours and secondary 
velocity profiles at the location of the vortex generator 
(x/D = 2.1). The location of the shed vortex along the 
downstream is shown in Fig 5.4. It is compared with the 
predicted location by the image vortex and the vortex 
interaction system. 

The total-pressure contours in Fig. 5.2 show that the 
boundary layer thickness in the region of downflow is 
decreased because the induced secondary flow pushes the 
high energy flow toward the wall. Adversely, the boundary 
layer thickness in the region of upflow is increased by the 
induced secondary flow. It shows that the appropriate 
vortex generators can control the main flow. 

The secondary velocity profiles in Fig. 5.3 show that 
the strength of the vortex decays in the downstream 
direction due to viscous diffusion. The streamwise vortex 
trajectory shows a helical character which is predicted by 
the inviscid theory. This characteristic is clearly shown 
in Fig. 5.4. 

Fig. 5.4(a) shows that the injected vortex moves 
radially inward. It shows significant deviation between the 
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computational result and simplified vortex interaction 
model. Even though the prediction model using the image 
vortex considers a mechanism which the injected vortex 
moves radially inward, it is very weak because the vortex 
moves radially inward after then the vortex core reaches 
the wall. The predicted locations obtained by two 
simplified prediction models are the same along the 
streamwise direction, as shown in Fig. 5.4(a). 

Physically, the boundary layer growth on the duct wall 
retards the growth of the vortex core to the wall, but the 
vortex core grows without blockage to the center of the 
duct. This causes a transverse pressure gradient which is 
not symmetric with respect to the vortex center. The 
pressure gradient between the vortex center and the duct 
wall is steeper than that between the vortex center and the 
center of the duct as shown in Fig. 5.2. The vortex moves 
radially inward as a result of this nonsynmetric pressure 
gradient. 

The location of the vortex along the circumferential 
direction is in agreement with the location predicted by 
the simplified model in the short region just downstream of 
the location of vortex generator as shown in Fig. 5.4(b). 
The deviation between the model location and computed 
results increases with increasing downstream distance. Even 
though the vortex interaction model considers the decay of 
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vortex strength by the wall effect, this model does not 
adequately consider the mechanism by which the vortex moves 
radially inward. At the same strength of vortex, if the 
vortex moves radially inward 10% from the original 
location, the induced velocity by the image vortex is 
reduced around 16% at the vortex core. The induced 
velocities overestimated by the simple models overpredict 
the azimuthal location of the vortex as it moves 
downstream. 

Fig. 5.5 shows the progression of the counter-rotating 
vortices of the same strength that rotate toward each other 
as they march down the duct. The left hand side and right 
hand side of Fig. 5.5 are the total-pressure contours and 
the secondary velocity profiles, respectively. Fig. 5.5(f) 
shows that the boundary layer thickness in the lower wall 
of the duct is one third of that in the upper wall of the 
duct at the station VI (x/D = 16.10). It shows that the 
main flow can be controlled by adjusting the number of 
vortex generators, strength and location of vortex, etc.,. 

The behavior of the vortices as they move downstream 
is qualitatively similar to the behavior predicted by the 
inviscid theory. Two vortices move away from each other, 
and also move radially inward as shown in Fig. 5.6(a). The 
deviation between the computational results and the 
predictions of the two simple models is due to the weak 
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mechanism of radial flow behavior in the simple models as 
mentioned in case of the single embedded vortex. 

The secondary velocities, induced by the counter- 
rotating vortices rotating toward each other, force the 
high energy flows into each other. Therefore, the pressure 
gradient with respect to the vortex axis in the case of 
counter-rotating vortices rotating toward each other is 
more symmetric than that of the single embedded vortex, 
whose larger induced velocities near the wall are 
associated with steeper pressure gradient near the wall. 
Comparing Figs. 5.4 and 5.6, one sees that the rate of 
radially inward motion when counter-rotating vortices are 
embedded as shown in Fig. 5.6 is approximately 7% lower 
than that when the single vortex is embedded as shown in 
Fig. 5.4. However, the rate of circumferential movement of 
counter-rotating vortices is lower than that of the single 
embedded vortex, even though the radial location of 
counter-rotating vortices is closer to the wall than that 
of the single embedded vortex. This is consistent because 
the secondary velocity induced by counter vortex acts 
oppositely to the direction which is induced on the vortex 
core by image vortex. 

Fig. 5.7 shows the progression of the counter-rotating 
vortices of the same strength that rotate away from each 
other as they march down the duct. The left hand side and 
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right hand side of Fig. 5.7 are the total-pressure contours 
and the secondary velocity profiles, respectively. Fig. 
5.7(f) shows that the boundary layer thickness at the lower 
wall of the duct is much greater than that at the upper 
wall of the duct. This is a contrary result compared with 
the case of the counter-rotating vortices of the same 
strength that rotate toward each other. 

The streamwise trajectories of vortices exhibit the 
same behavior as that predicted by the inviscid theory. The 
vortices attract each other in a short region downstream of 
after the vortex generators, and then they proceed to march 
away from the wall. As the two vortices move closer to each 
other, the pressure gradient between the vortex center and 
the duct wall is increased, but the pressure gradient 
between the vortex center and the center of the duct is 
decreased. Downstream of station IV (x/D = 8.20), the 
pressure gradient between the vortex center and the 
symmetric line of two vortices is steeper than that between 
the vortex center and the duct wall as shown in Fig. 5.8. 
It means that the two vortices move away from each other 
during the time they proceed to march away from the wall as 
shown in Fig. 5.7. 

Fig. 5.7 shows that the predicted vortex location by 
two models is overpredicted except in a short region 
downstream of the vortex generators. This deviation occurs 
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from the weak mechanism of the two simplified models as 
previous mentioned. Two prediction models do not have a 
mechanism which each vortex tries to settle at a stable 
location, i.e., vortex moves to the position of radially 
symmetric pressure gradient. 

Fig. 5.10 shows the total-pressure contours along the 
duct when the co-rotating vortices are embedded; the 
secondary velocity profiles are shown in Fig 5.11. The 
strength of vortex (A) is twice the strength of vortex (B) . 
As the vortices march down the duct, the circumferential 
movement of vortex (A) is faster than that of vortex (B) 
because of its large induced velocity on the vortex core. 
This is anticipated by the inviscid theory. The vortex (B) 
is collapsed into the vortex (A) at some distance as shown 
in Fig. 5.11 because two vortices have the same direction 
of vorticity. 

Figs. 5.12 and 5.13 are the locations of vortex (A) and 
vortex (B) along the duct, respectively. They agree well 
with the results by the prediction models in a short region 
just after vortex generators. The deviation between the 
computational results and the prediction by two models 
occurs from the weak mechanism of the two prediction models 
as previously mentioned. 

The satisfactory results of the computation in the 
straight duct with vortex generators suggest that the 
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vortex model employed in this work can be applied to solve 
the full three-dimensional Navier-Stokes equations. The 
internal flow can be controlled by an appropriate 
adjustment of the location, strength, lateral spacing, and 
number of vortex generators. The computational results 
agree well with the results of the prediction models in a 
short region just after the location of the vortex 
generators, even though we adopted the same constants which 
were derived from the experimental results on the external 
flow. For the better prediction of vortices along 
downstream in the internal flow, an experiment in the duct 
with vortex generators is necessary to find the correct 


constants . 





Fxg. 5.1(a) Computational grid for the vortex interact 
studies within a cylindrical duct, L/D = 20 
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5.1(b) Measurement stations along the circular 
straight duct. 
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(d) station IV (x/D = 8.20) 


110 



(f) station VI (x/D = 16 . 20 ) 











(d) station IV (x/D = 8.20) 




Prediction by two models 
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prediction by vortex interaction system 
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Normalized distance along duct (.\/D) 

Fig. 5.4(a) Radial trajectory of the single embedded 
vortex. 
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Fig. 5.4(b) Angular trajectory of the single embedded 
vortex. 
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Fig. 5.6(a) Radial trajectory of the counter rotating 
vortices of the same strength that rotate 
toward each other. 
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Fig. 5.6(b) Angular trajectory of the counter rotating 
vortices of the same strength that rotate 
toward each other. 
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STATIC PRESSURE COEFFICIENT CONTOURS 
( Upflow pairs ) 




(b) station IV (x/D = 8.20) 

Fig. 5.8 Static-pressure coefficient contours of the 

counter-rotating vortices of the same strength 
that rotate away from each other. 

( Ma = 0.6/ Re d = l.OxlO 4 ) 
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Fig. 5.9(a) Radial trajectory of the counter-rotating 

vortices of the same strength that rotate away 
from each other. 
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5.9 (i) Angular trajectory of the counter-rotating 

vortices of the same strength that rotate away 
from each other. 
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TOTAL PRESSURE COEFFICIENT CONTOURS 
{ Co-rotating vortices ) 
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Fig. 5.10 Total -pressure coefficient contours when the 

co-rotating vortices are embedded. Vortex (A) has 
double the strength of vortex (B) . 

( Ma = 0.6, Re d = 1.0X10 4 ) 


(c) station III (x/D = 4.24) 



(d) station rv (x/D = 8.20) 
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(£) station VI (x/D = 16.20) 




(b) station II (x/D = 2.90) 


Pig. 5.11 Secondary velocity profiles vhen the 

co-rotating vortices are embedded. Vortex (A) has 
double the strength of vortex (B) . 

( Ma = 0.6, Re. = 1-OxlO 6 ) 







(e) station V (x/D = 12.20) 



(f) station VI (x/D = 16.20) 
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Fig. 5.12(a) Radial trajectory of the vortex(A) which has 
double the strength of the vortex (B) in the 
co-rotating vortices. 
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Fig. 5.12 (b) Angular trajectory of the vortex (A) which has 
double the strength of the vortex (B) in the 
co— rotating vortices - 
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Fig. 5. 13 (a) Radial trajectory of the vortex(B) which has 
half the strength of the vortex (A) in the 
co-rotating vortices. 
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7ig. 5. 13 (Is) Angular trajectory of the vortex(B) which has 
half the strength of the vortex (A) in the 
co-rotating vortices. 



CHAPTER 6 


DIFFUSING S-DUCT WITH VORTEX GENERATORS 


6.1 Vortex Model 

The shed vortex from the vortex generators is modeled 
by introducing two-dimensional vortical flow in the cross- 
plane. Eq. (5.10) provides this vortical structure. 


u 9 = 


U I 2 

[ 1 - exp ( - — ^ — ) ] 


(5.10) 


r is the vortex strength at the tip of the vortex 
generator. The r term is a function of the geometry of the 
generator, and the oncoming flow conditions, r is defined 
by the strength of one fully rolled up trailing vortex; 

T = C. c, u (6.1) 

• 2 1 


where c L is the lift coefficient, c, is the chord length of 
the vortex generator, u is the velocity of the flow at the 
generator tip, and C a is the constant which considers the 
viscosity and turbulence effect, etc.,. C a cannot be greater 
than 0.45 according to inviscid wing theory and by 
experiment. 
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Three pairs of one half of the NACA0012 wing section 
were used as the vortex generators in the Vakili et 
al.(1986) experiment. They were installed in the duct at 
S/D; = 0.087, and at azimuthal angles of -41.4 s , 0.0° and 
+41.4. The height and chord length of the vortex generator 
were h/D t = 0.715X10' 1 and q/Dj = 0.103, respectively. The 
vortex generator pairs had geometric incidence angles of 
+14° and -14° relative to the duct centerplane. 

Eq. (6.1) can be expressed in nondimensional form; 


r = c — — — 

Dil L ' “ 2 D i U m 


( 6 . 2 ) 


From the experimental conditions, u/U. is taken as 1 and c L 
is assigned as 1.4 because the incident angle of the vortex 
generator is 14°. In this study, six different vortex 
strengths r/DjU. = 0.005, 0.010, 0.015, 0.020, 0.025 and 

0.03 0 are investigated to compare with the experimental 
data and to study the parametric effect of different vortex 
strengths. When the vortex strength (r/DiU.) is equal to 
0.030, C a is 0.4. In the choice of various vortex strength, 
the decreasing of the vortex strength implies that the 
incident angle of the mounted vortex generator is 
decreasing. 

In Eq. (5.10), the length x is estimated to be the 
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distance 0.087D; from the location of the vortex generator 
to the crossplane of the vortical structure. The rotational 
velocities at the cross plane (S/D, = 0.17) are evaluated 
with the image vortices based on the circle theorem 
mentioned in section 4.2.1. If the rotational velocity near 
the vortex core is greater than U^/5, the velocity at that 
point is assumed to be of that magnitude. These secondary 
velocities of the vortex model are added to the secondary 
velocities without vortex generator at the same plane. The 
combined vortical structures are applied as the source term 
in the crossplane. Fig. 6.1 shows the secondary flow 
structure at this plane. 

In this computation, the residuals for these numerical 
solutions were reduced approximately three orders of 
magnitude. Solutions were obtained on the Cray-YMP. The 
number of iterations required to obtain the converged 
solutions were approximately 25,000. The computational 
speed was 950 iterations per CPU hour. The mass flow 
changes between the inlet and exit was within 1 percent. 


6.2 Results and Discussion 

Fig. 6.2 shows values of the computed static-pressure 
coefficients (continuous curves in the figure) and 
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experimental values (symbols) at — 10° , 90° and 170° for 
various vortex strengths. Numerical results with T/D.U,, = 
0.025 are close to the experimental data. In the following 
discussion, these numerical results will be compared with 
the experimental results obtained by Vakili et al.(1986). 
The surface static-pressure distribution at <p - 170° shows 
some deviation between the experimental and numerical 
results near the location of the vortex generators. Recent 
experimental results on the same geometry by Reichert and 
Wendt (1992) show that there is no perceptible upstream 
influence on the static-pressure distribution, caused by 
the vortex generator arrays. In their experiment, Wheeler 
wishbone generators are used. This type of generator forms 
a pair of counter-rotating vortices with the flow between 
vortices directed upwards. However, the experimental data 
obtained by Vakili et al. (1986) show some influence on the 
static-pressure distribution by the vortex generator 
arrays. The influence of the vortex generator arrays on the 
static-pressure distribution depends on the location of the 
vortex generators and data acquisition points, but the 
vortex model employed in this study shows very little 
upstream influence on the static-pressure distribution, as 
can be seen in Fig. 6.3. 

For the smallest vortex strength (T/DjU. = 0.005) , the 
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results differ only marginally from the flow without vortex 
generators as shown in Fig. 6.3(a). Figs. 6.3(b) and 6.3(c) 
show higher values of static-pressure for the larger vortex 
strength in the second half bend of the duct. The static- 
pressure distribution lines cross each other at the 
inflection point of the duct (S/D; = 2.62). The static- 
pressure value at the cross point is less than the peak 
value at <p = 0° near S/D; = 2.5. These results are very 
similar to the experimental results conducted by Reichert 
and Wendt (1992) . 

In Figs. 6.3(a) and 6.3(b), the constant static- 
pressure values at <p = 170° in 3 < S/D; < 4 are associated 
with the flow separation. These figures show that the 
region of the constant static-pressure value decreases with 
increasing the vortex strength. The reverse flow of 
streamwise velocity dose not occur when the injected vortex 
strength is greater than r/DjU,, = 0.020. 

Fig. 6.4 shows the secondary velocity profiles 
compared with the experimental results. The computed and 
experimental results show on the right hand side and left 
hand side, respectively; only half of the cross-plane is 
shown because the flow is symmetric along the duct cross 
section. The numerical results agree closely with the 
experimental results except the behavior of the vortices (C) 
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along the downstream. Note that the resolution and the 
locations of data collection points of the experiment and 
computations are different. The computational results of 
Fig. 6.4(a) are plotted in a denser resolution in order to 
show more clearly the vortices just downstream of the 
vortex generators. 

At the first half bend, the high energy flow moves 
toward the upper wall and the low energy flow migrates 
circumferentially from the upper wall to the lower wall. In 
Fig. 6.1, the rotational velocities of the injected 
vortices (B) have the same direction as the low energy flow 
near the wall, but vortices (A) and (C) have opposite 
rotational velocities to the low energy flow. This makes 
the secondary velocities of vortices (B) results in stronger 
than those of the other vortices. It also makes vortices (C) 
more quickly decaying. In the experimental results, the 
vortices (C) do not decay as quickly as in the computation; 
even if the strength of vortices (C) is weaker than the 
other vortices. The low energy flow at the vortex plane is 
retarded by the installed vortex generators on the wall. 
This means the injected vortices have little influence from 
the low energy flow. The vortical structure of the vortex 
model is strongly influenced by the low energy flow at the 
location of the vortex generators. 

Fig. 6.5 shows the total-pressure contours compared 
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with the experimental results. The small effect of the 
vortices (C) is clearly shown in this figure. The thickness 
of the computed "boundary layer" at <p = 90° is less than 
observed experimentally. 

The static-pressure contours are shown in Fig 6 . 6 . The 
numerical results agree qualitatively with the experimental 
results. In the first half bend, higher static pressure is 
shown near the upper wall because of the duct curvature. 
Opposite behavior is shown in the second half bend owing to 
the same reason. 

The variation of the boundary layer thickness at 93 = 
10°, 90° and 170° along the duct is shown in Fig. 6.7. The 
boundary layer thickness is defined as the normal distance 
from the wall where the total-pressure coefficient is 1 . 0 . 
The boundary layer thickness of the flow with vortex 
generators depends greatly on the vortex strength. The 
computed boundary layer thickness at <p - 90° is less than 
the experimental result because the injected vortices (C) 
are weaker than the experimental values. However, the 
computed results show that the trend of the boundary layer 
thickness variation along the duct is quite similar to the 
experimental results. 

Fig. 6.8 shows the total-pressure contours with and 
without vortex generators. The right hand side and left 
hand side show the numerical results with and without 
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vortex generators, respectively. The effect of the injected 
vortices is clearly shown in this figure. The injected 
vortices push the high energy flow toward the lower energy 
region. This resulting force prevents the flow separation 
at the inflection point of the duct. 

Fig. 6.8(a) to (c) show that the boundary layer 
thickness of the flow with vortex generators near the upper 
wall is less than that of the flow without vortex 
generators. This results from satisfying a constant mass 
flux because the shed vortices injected near the lower wall 
cause a streamwise velocity deficit in the region of the 
vortex core. However, the experimental results do not show 
any difference between the boundary layer thicknesses with 
and without vortex generators in the upper wall of the 
first half of the duct. This probably results from 
deficiencies in the experiments, primarily coarse data 
acquisition locations and uncertainties in the total 
pressure measurements using pitot tubes. 

Fig. 6.9 shows the secondary velocity profiles with 
and without vortex generators. The interaction between the 
injected vortices and the counter-rotating vortices 
resulting from the flow separation is clearly shown in this 
figure. The injected vortices suppress the growth of these 
counter-rotating vortices. 

The static-pressure contours with and without vortex 
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generators are shown in Fig. 6.10. At the vortical plane 
(S/D, = 0.17), the distortion of the constant static- 
pressure contours is a result of the injected vortices. The 
change of the static-pressure along the duct shows the same 
flow phenomena as mentioned in the discussion of the flow 
without vortex generators. Figs. 6.10(d) and 6.10(e) show 
that the constant static-pressure contours are flatter in 
the low energy flow region of the second half bend. The 
injected vortices result in a more uniform flow and higher 
diffusion at the exit than occurs without vortices. 

Figs. 6.11 and 6.12 show the numerical results with 
the vortex strength (T/DiU* = 0.015). The computed results 
show that the effect of the injected vortices is weaker 
than with the strong vortex strength (T/DiU^ = 0.025), as 
one could expect. The region of diminished total-pressure 
at the exit is larger and the static-pressure contours are 
more distorted. Fig. 6.12 shows the interaction between the 
injected vortices and the counter-rotating vortices 
resulting from the flow separation. It also shows that the 
growth of these counter-rotating vortices are suppressed by 
the injected vortices. Fig. 6.12(e) shows that the 
secondary velocities between the vortices (A) are 
overestimated. This results from the small eddy viscosity 
in the flow separation region as mentioned in the 
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discussion of the flow without vortex generator. 

Fig. 6.13 shows the total-pressure contours at the 
exit with different vortex strengths. The region of 
diminished total-pressure is significantly reduced with 
increasing the vortex strength. 

A total-pressure recovery (C,J is calculated using 
area weighted values from the computational mesh over the 
cross stream plane. 


Cp. ’ J / C pc dA 


Using a similar method, the total-pressure recovery of a 
segment is determined by integrating the total-pressure 
coefficient over a segment of the cross stream plane of 
angular extent <p. 


C po (cp) 


fdA 


( 6 . 4 ) 


A distortion coefficient is useful to describe the 
efficiency of inlet duct or to compare the performance of 
several inlet ducts. There are many ways to define the 
distortion coefficient depending on the comparison 
purposes. Distortion coefficients measuring radial or 
circumferential distortion have been used. Early workers 
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simply defined the distortion coefficient in experiment to 
be the difference of the normalized maximum rake total- 
pressure and the normalized minimum rake total-pressure. In 
this study, the distortion coefficient DC(^) is defined 
using the cross stream plane segment that results in the 
lowest value of . In case of S-duct, the segment angle 
ip is defined to the azimuthal angle from the centerline in 
the lower energy region. The values of p are chosen to 60°, 
90° and 120°. 

DC ( <p ) = C po - C po ( <p ) (6.5) 

Fig. 6.14 show the total-pressure recovery at the exit 
with various vortex strengths. Fig. 6.15 show the 
distortion coefficient at the exit. For the smallest vortex 
strength (T/DjU^, = 0.005), the total-pressure recovery is 
slightly reduced. This indicates that the small vortex 
strength acts as flow blockage at the location of the 
vortex generators. The vortex strength is quickly reduced 
in the first half bend. The resulting force is not enough 
to suppress the counter-rotating vortices resulting from 
the flow separation. Small vortex strength is seen to 
affect the flow in a detrimental way. This phenomenon with 
small vortex strength is shown in the experiment by 
Reichert and Wendt(l992). The total-pressure recovery and 
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distortion coefficient significantly improved 
increasing vortex strength. 
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SECONDARY VELOCITY VECTORS 
( Vortex strength, r/DiU. = 0.025 ) 
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SECONDARY VELOCITY VECTORS 


With Vortex 
Generators 
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Pig. 6.9 Comparison of the secondary velocity profiles 
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CHAPTER 7 


CONCLUSION AND RECOMMENDATION 


The numerical results on a diffusing S-duct without 
vortex generators show the phenomena of three-dimensional 
flow separation. The computed results agree well with the 
experimental results except in the flow separation region. 
Downstream of flow separation, the strength of the 
streamwise velocity deficit (U a - u) is overestimated at 
the region of the counter-rotating vortices resulting from 
the flow separation. This results from underestimating the 
eddy viscosity effect in the flow separation region by the 
turbulence model. However, the computed results are better 
than the previously published work obtained by Harloff et 
al. (1992b) with an alternative turbulence model. In order 
to obtain better solutions in the flow separation region, 
further efforts on three-dimensional turbulence modeling 
are necessary. 

The computed results on a straight duct with vortex 
generators show how the injected vortices decay, move along 
the duct, and interact with the boundary layer in a simple 
geometry. For a short region (approximately three times 
diameter) downstream of the vortex generators, the vortex 
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core locations determined from the simplified model and by 
the full computations are in good agreement. Farther 
downstream, however, the simplified model is not able to 
predict the radially inward motion. In order to provide 
more accurate vortical structure for the vortex generator 
model, experiments with vortex generators in straight ducts 
would be useful. 

The computed results on a diffusing S-duct with vortex 
generators show the interaction between the separated flow 
and the injected vortices. As the strength of the vortex 
generators increases, the extent of flow separation region 
is decreased. When the strength of the injected vortex is 
greater than r/D.U^ = 0.020, reverse flow along the 
streamwise direction does not occur. 

The computed results depict well the behavior of the 
injected vortices as they travel downstream except for the 
injected vortices that are introduced into the region with 
strong secondary velocity induced by the curvature of the 
duct. The behavior of the injected vortices along the 
streamwise direction depends on the induced secondary 
velocity and the injected location within the duct, even if 
the vortices are injected with the same strength. 
Experiments are needed to obtain the secondary velocity 
just downstream of the vortex generators in order to obtain 




180 


an accurate vortical structure for modeling the shed 
vortices in the curved duct. 

The total-pressure recovery increases and the 
distortion coefficients decrease at the exit with 
increasing vortex strength, except for the smallest vortex 
strength (T/DjU^ = 0.005). This indicates that there exists 
an optimal vortex strength which will minimize the flow 
distortion at the exit. In order to obtain the optimum flow 
at the exit, additional numerical studies are necessary 
with various axial locations, lateral spacing, height, and 
number of vortex generators. 
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APPENDIX A 


The Elements of Jacobian Matrices 
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where g ~ ( u2 + v 2 + w 2 ) / 2 

The elements of B and C are similar to those of A, 
only | in the element of A has to replace to tj and f, 
respectively. 



APPENDIX B 


Vortex Trajectory In a Tube Using the Image Vortex System 


In order to form a simple model which estimates the 
trajectory of a vortex in a tube, the two-dimensional 
problem of the motion of an infinite line vortex in a 
circular cylinder is considered, and then superpose an 
axial velocity to describe the motion of a traveling vortex 
in a circular cylinder. 

Consider first the inviscid model. Fig. B.l is a 
diagram of the vortex in a tube with an image vortex. X = 
r/R is the non-dimensional radial location of the vortex. 
Also note that the image flow must have a superposed 
circulation. This plays no role in the following 
discussion. The internal vortex moves with the velocity 
induced by the image vortex: 


T X 
2^ 1 -X 2 


(B.l) 


The angular velocity of the vortex motion is u,/\R. Thus: 


(0 


r i 

2* R 2 1 - X 2 


(B.2) 


Now consider a decaying vortex. The analytic solution for 
a decaying vortex at the origin, with laminar flow is: 
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u a = 


r -Ur 2 

L [ 1 - exp ( - ) ] 
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4 v x 


(B • 3 ) 


For the axially-moving vortex, the time is the particle 
travel time for axial motion: 


t = 


x 

U m 


(B. 4) 


The azimuthal velocity is then given by: 


u e = 


2 tc r 


[ 1 - exp ( 


4 v r 


) ] 


( B . 5 ) 


If identically decaying vortices are employed in the image 
system, the tube boundary is no longer a streamline. 
However, the image system remains valid using a simple 
approximation to the flow field described by Eq. (B.5). The 
decaying vortex has a finite velocity slope at r = 0, and 
behaves like an inviscid vortex as r — ». it can therefore 
be approximated as an inviscid vortex with a solid body 
core. The radius of the core is the value where the solid 
body velocity and the vortex velocity are equal. The core 
radius for a vortex at the origin is: 
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For a vortex whose center is at \R, the outer radial 
location of the solid body core is given by: 


r 


IT? 


X R + 2 


Z* = R 


\ 

x/R 

u. 

1 'I 

U m R / v J 


( B . 7 ) 


The core reaches the wall of the tube where r m = R. This 
occurs at an axial location given by: 



UR 


(1 - X) 


(B. 8) 


For x/R < ( x/R) a , the azimuthal velocity and angular 
velocity of the vortex are given by the first factors in 
Eq . (B. 5) : 


u e = 


r x 

2-x R {i - X 2 ) 


(B . 9 ) 


Q 


r i 

2- R 2 (1 - X 2 ) 


(B. 10) 


To this point the model gives no information about the 
radial migration of the vortex. We now use the growth of 
the solid body core to obtain an estimate, albeit weak, of 
the migration of the vortex center toward the tube center. 
If we think of the core as a solid body, then continued 
growth of the core beyond x = x a forces migration of the 
vortex toward the tube center. Then Eq. (B.8) gives the 
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radial location of the center of the vortex for any x > x a , 
and Eqs. (B.9) and (B.10) give the corresponding azimuthal 
and angular velocities. A convenient measure is the angle 
of travel of the vortex: 


- f cj dc ~ ^ ( 

r *. dx 

* f x dx 'l 

j 2 n R 2 U m l 

1 - xl 

J*. 1 - l 2 , 


(B.ll) 


6 = 


2 n R 2 n 


l - K 


r * dx ' 
■U l - X 2 


( B . 12 ) 


In the integral in Eq . (B.12), if Eq. (B.8) holds without 

the subscript "a", one then obtains: 


1-1 


N 


4 xj_R 
U m R/ v 


s /c ? 


(B. 13) 


Eq. (B.12) then becomes 


r 

1 U m R 2 

(° do ' 

2 it R 2 U m 

i - J 2 4 v J 

°* Jo ( 2 - /o) y 


(B . 14 ) 


Performing the integration yields: 


r 

x - . 1 i n 

/ 

1 



2 tz R z U m 

1 - A 2 . 2 v 

1 

- -I/O 




V 

2 J 



But 
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and 

Therefore : 


Eqs . (B . 13) , 

For a < ( 1 - 
is given by: 


/°7 = 1 - 


(B.16) 


i u m R 2 

4 v 


(1 - X,) 2 


(B . 17 ) 


0 = 


8itv 


1 - *c 

1 + X. 


- 2 In 


2 ~ y/~Q 

1 + X. 


(B. 18) 


X = 1 - fa 


f5 * 



xJ_R 
(J m R / v 


(B . 13a) 
(B. 13b) 


B.13a) and (B.13b) hold in the range: 

( 1 - X a ) 2 s o s 1 (B . 19 ) 

y. a ) : , X = X a and the azimuthal angular travel 


0 = 


Tx 

2 n R z U m ( 1 - X 2 ) 


(B.20) 


An alternative measure of the azimuthal travel is the 
tangent of the helix angle. This is obtained by merely re- 
writing Eqs. (B.18) and (B.20): 
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OR 

X 


r \ 

2 t i RU m o 


1 ~ *a 

1 + 


( 

- 2 In 

k 


2 - Vo ) 
1 * A J 


(B . 2 1 ) 


b_r = r 

x 2 it R U m ( 1 - A. 2 ) 


(B.22) 


where Eqs . (B.21) and (B.22) are valid in the large a and 

small cr ranges, respectively. 
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Fig. B.l When a single vortex of strength r is located 
inside the circular cylinder, image vortex is 
located on the line connecting the center of a 
circular cylinder and a single vortex. 


APPENDIX C 


Vortex Trajectory In a Tube Using the Vortex Interaction 

Model 

Circulation decay of vortex in the turbulent flow is 
faster than in the laminar flow because of a large eddy 
viscosity. Circulation decay on the previous model is very 
slow because only kinematic viscosity is used. In this 
model, the wall effects and proximity effects are 
considered to predict circulation decay of the vortex in 
the tube. At first, circulation decay by the wall effects 
is considered. Fig. C.l is a diagram when a vortex is 
embedded at some crossplane location x. The secondary flows 
produced by the injected vortex give rise to a 
corresponding circumferential component of the wall shear 
stress (r ri ). In turn, this stress results in a torque 
opposing the rotation of the vortex. The moment M i opposing 
the rotation of vortex i can be obtain by integrating the 
magnitude of the elemental torque: 

M i = a | f 1 c© dx , f ) (C.l) 

where |r| is the distance from the center of vortex to the 
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wall, and H t is the angular momentum of vortex i of 
elemental thickness dx . 

Assuming proportionality between the vortex angular 
momentum and its circulation, we obtain the streamwise 
circulation gradient: 


dT 

dx 


1 

CW P 4. 



f 1 d 0 , 


f ) 


( C . 2 ) 


where C l(/ is the unknown constant of proportionality with 
units of m 2 . 

is a function of wall coordinates, i.e., r rt = 
r rt (r, 0) and it is proportional to the circumferential 
component of the secondary velocity at the wall. To 
simplify this expression, the correlation suggested by 
Pauley and Eaton (1988) is adopted. 

*re = C t *zd u * ( 0 ' r = | f | ) (C.3) 


where u 0 (6,r = |r|) is the circumferential velocity with 
image vortex at the wall, r,* is the wall shear stress of 
the corresponding two dimensional boundary layers and C, is 
a scaling factor with units of sec/m; 


t 


2D 



P m uj 

2 


_ JL 

( 0.3164 Rs d 4 ) 


(C.4) 


Re d is a Reynolds number based on the duct diameter. 
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Substitute Eqs. (C.3) and (C.4) into (C.2), and then we 
obtain the circulation gradient along the streamwise 
direction by the wall effects; 


dT 

dx 


c, u m 
~c~7 


( 0.15 82 Re d 4 ) /* 2,t |f| 2 cos(p)c?9 

J 0 


(C. 5) 


where 0 is an intersectional angle between the r-direction 
and the 5-direction. 

When the multiple vortices are embedded in a 
crossplane, the vortices interact both because of their 
induced field and through diffusion. Fig. C.2 is a diagram 
of two counter-rotating vortex cores in close proximity. 
The circulation decay by proximity effect can be expressed: 


dT | — o ^ | do I 

dx 1 prox p:ox | r | 1 dR lr * r ° 


(C.6) 


where C prac is the unknown constant of proportionality with 
units of m 2 , and the sign of C prai depends on the rotation 
direction of the neighboring vortex. 

Total circulation decay is written: 


dT = dT | + dT | 

dx ~ dx [wi dx lprox 


(C.7) 


With this gradient and an assumption that the embedded 
vortices move axially at the free streamwise velocity, we 
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can track the streamwise displacement of the embedded 
vortex step by step 

r (x+ix) = r ( x ) + ix ( — } — ) (c. 8) 

c X 

0 ( x ->• A x ) = 0 ( x ) » Ax ( - 9 . { ) (C.9) 

OX 

r(x + Ax) =r(x) i Ax( * n_x) ) (c.io) 

a x 

As noted early, the same constants C , (1. 40x10^ m 2 ) , 
C prat (1.40X10 -4 Til- ) and C T (0.046 sec/m), which were derived 
from the experimental result on the flat plate, are adopted 
to predict the circulation decay in the internal flow. 




Fig. C.l The secondary flow field generated by vortex i 
gives rise to a local circumferential component 
of wall shear stress which opposes the 

rotation of vortex i. 



Fig. C.2 Two neighbor vortices in a circular duct for 

evaluating proximity circulation losses. The r 
represents the coordinate axis along the line 
connecting adjacent cores, and r 0 is the location 
on r where the vorticity changes sign in the 
model . 
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